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It  is  a  well-established  fact  that  the  state  of  the  ground  is  driven  in  a  large  part  by  the  downwelling 
solar  and  IR  fluxes.  Models  developed  to  predict  the  state  of  the  ground  depend  critically  on  these  fluxes 
for  initialization.  When  measured  solar  and  infrared  fluxes  are  not  available  they  must  be  computed.  We 
have  compared  the  ground  temperatures  as  computed  by  the  thermal  model  SWOETHERM  using  differ¬ 
ent  solar  flux  initialization  schemes.  These  initialization  schemes  used  measured  solar  flux  values  ob¬ 
tained  during  the  Smart  Weapons  Operability  Enhancement  (SWOE)  field  programs,  and  calculated  solar 
flux  values  from  a  semi-empirical  model  (Shapiro’s  model),  a  plane  parallel  model  (MODTRAN),  and 
ARL’s  AIM  (Atmospheric  Illumination  Module)  model.  We  investigated  the  response  of  the  surface  tem¬ 
perature  to  different  solar  flux  initialization  schemes  while  all  other  environmental  parameters  were  held 
constant.  We  found  that  for  clear  skies  all  schemes  resulted  in  nearly  identical  surface  temperatures.  For 
partly  cloudy  and  cloudy  skies  only  the  AIM  model  can  mimic  the  spatial  variability  observed  with  the 
measured  solar  fluxes.  The  Cloud  Scene  Simulation  Model  (CSSM)  was  used  to  determine  the  spatial 
variability  of  the  clouds.  The  cloud  distributions  were  then  used  by  AIM  to  produce  the  variations  of  the 
surface  solar  loading.  CSSM  also  has  the  capability  to  produce  the  temporal  variations  in  the  cloud  fields 
for  short  periods  of  time.  Thus,  it  would  be  possible  to  use  CSSM  and  AIM  to  produce  the  temporal  and 
spatial  variations  in  the  solar  loading.  Models  like  AIM  frequently  incur  a  large  computational  burden.  In 
order  to  reduce  the  computational  burden  associated  with  AIM  we  have  implemented  several  new  proce¬ 
dures.  Distributed  energy  budget  models  used  to  predict  the  state  of  the  ground  require  distributed  envi¬ 
ronmental  parameters  for  initialization.  Many  of  these  parameters  can  be  obtained  from  mesoscale  mod¬ 
els  such  as  MM5  or  databases  associated  with  programs  such  as  IMETS.  But,  to  our  knowledge,  none  of 
these  models  or  programs  provide  distributed  solar  or  infrared  fluxes,  a  key  initialization  parameter  of 
energy  budget  models.  Models  like  AIM  linked  to  CSSM,  or  for  that  matter  any  model  that  provides 
information  on  the  spatial  and  temporal  distribution  of  atmospheric  conditions,  can  be  used  to  provide  the 
spatial  and  temporal  distribution  of  radiative  fluxes. 
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1  INTRODUCTION 

It  is  a  well-established  fact  that  the  state  of  the  ground  is  driven  in  large  part 
by  downwelling  solar  and  IR  fluxes.  Models  developed  to  predict  the  state  of  the 
ground  for  Army  operations  will  depend  critically  on  these  fluxes  for  initializa¬ 
tion.  Unfortunately,  these  fluxes  are  not  routinely  measured  parameters  as  is  the 
case  with  more  common  meteorological  parameters  such  as  temperature  and 
relative  humidity,  so  indirect  methods  must  be  used  to  generate  the  required  flux 
initialization  information  for  state-of-the-ground  models.  One-dimensional 
models  used  in  a  distributed  mode  to  predict  the  state  of  the  ground  over  some 
spatial  domain  require  distributed  initialization  parameters  over  those  domains. 
Assuming  the  atmosphere  overlying  the  spatial  domain  of  interest  is  horizontally 
homogeneous  may  lead  to  erroneous  results.  This  would  be  especially  true  of 
partly  cloud  sky  conditions  where  the  solar  loading  can  vary  significantly  over 
fairly  small  spatial  extents.  Davis  and  his  colleagues  (1998)  have  developed 
surface  segmentation  techniques  based  on  terrain  topography  (slope  and  aspect) 
and  land  surface  characteristics  (soil  type,  vegetation  cover,  etc.).  They  found 
that  segmenting  a  1  -  by  1  -km  area  at  Grayling,  Michigan,  into  250  regions  and 
running  a  one-dimensional  thermal  model  for  each  region  enabled  them  to 
accurately  reproduce  the  spatial  and  temporal  distribution  of  the  snow  cover. 
Mesoscale  models  such  as  MM5  have  better  forecast  skills  at  12  and  36  km  than 
at  higher  resolutions.  A  resolution  of  12  km  may  be  insufficient  for  the  initiali¬ 
zation  of  distributed  state-of-the-ground  models,  especially  for  radiative  flux 
parameters  and  segmentations  on  the  scale  that  Davis  found  to  be  optimal.  This 
raises  the  question  of  the  importance  of  using  distributed  initialization  parameters 
for  distributed  modeling  of  the  state  of  the  ground.  We  will  investigate  the  sen¬ 
sitivity  of  the  ground  temperature,  a  state-of-the-ground  parameter,  as  predicted 
by  the  S  WOETHERM  thermal  model  to  different  flux  initialization  techniques, 
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including  initialization  using  measured  flux  values  obtained  during  the  Smart 
Weapons  Operability  Enhancement  (SWOE)  field  programs. 

The  Cold  Regions  Research  and  Engineering  Laboratory  (CRREL)  has  nu¬ 
merous  data  sets  that  can  be  used  to  initialize  the  solar  flux  models.  One  of  the 
more  comprehensive  ones  was  collected  during  the  JT&E  SWOE  program.  These 
data  sets  also  contain  the  information  that  can  be  used  as  ground  truth  for  the 
evaluation  of  predicted  solar  fluxes.  Solar  fluxes  will  be  calculated  using  a  semi- 
empirical  scheme  developed  at  CRREL  based  on  the  work  of  Shapiro  (1972, 

1 982, 1 987),  a  plane-parallel  scheme  using  MODTRAN,  and  ARL’s  AIM 
(Atmospheric  Illumination  Module)  model.  Only  the  AIM  model  produces 
distributed  flux  values  over  a  user-defined  spatial  domain.  AIM  uses  the  Cloud 
Scene  Simulation  Model  (CSSM)  (Cianciolo  and  Rasmussen  1996)  in  conjunc¬ 
tion  with  the  Boundary  Layer  Illumination  and  Transmission  Simulation  (BLITS) 
radiative  transfer  program  to  determine  the  spectral  and  spatial  distribution  of 
fluxes  in  cloudy  and  clear  atmospheres.  Unlike  the  first  two  approaches  that  are 
either  a  parameterization  or  assume  a  plane-parallel  atmosphere,  BLITS  uses  a 
physics-based  approach  that  models  three-dimensional  fluxes  through  dense 
clouds.  Scenarios  will  be  run  using  data  sets  from  Yuma  and  Grayling  for 
different  sky  conditions  (clear,  partly  cloudy,  and  overcast). 
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2  SOLAR  RADIATION  INITIALIZATION  METHODS 


AIM  surface  radiation  calculations 

The  Atmospheric  Illumination  Module  has  been  extended  to  permit  surface 
solar  loading  calculations  in  addition  to  basic  surface  illumination  calculations. 
AIM  was  originally  an  extension  of  previous  work  accomplished  as  part  of  the 
ARL  WAVES  (Weather  and  Atmospheric  Visualization  Effects  for  Simulations) 
program  to  permit  better  computation  of  propagation  and  illumination  effects  in 
the  presence  of  dense  clouds.  The  AIM  propagation  model  BLITS  leveraged  off 
useful  properties  of  the  WAVES  model  BLIRB  (Boundary  Layer  Illumination 
and  Radiation  Balance),  using  the  discrete  ordinates  method  (DOM)  for  charac¬ 
terizing  angular  integration  and  scattering  effects  due  to  atmospheric  haze  aero¬ 
sols.  However,  BLIRB  did  not  perform  well  when  dense  aerosols  (cumulus  or 
stratus-type  clouds)  were  present  in  the  modeling  volume.  Thus  BLIRB  could  not 
adequately  characterize  partly  cloudy  conditions.  The  BLIRB  model  had  been 
designed  to  divide  the  modeled  volume  (typically  a  volume  several  kilometers 
deep  starting  at  the  earth’s  surface)  into  a  series  of  rectilinear  cells.  In  three 
dimensions  (3-D),  these  cells  were  then  populated  with  haze  and  cloud  aerosols. 
BLITS  uses  this  same  technique,  plus  an  improved  technique  for  representing  the 
scattering  properties  of  cloud  aerosols:  EOSAEL’s  CLTRAN  (CLoud  TRANs- 
mission)  Khrgian-Mazin  model  for  cloud  aerosols  and  its  predictions  of  vertical 
structure  of  aerosol  density  and  properties  within  a  cloud  layer.  In  addition,  a  log- 
least-squares  technique  was  adopted  to  determine  the  Legendre  expansion  coeffi¬ 
cient  used  for  representing  the  cloud  aerosols  input  parameters  to  the  radiative 
transfer  model  BLITS.  Due  to  the  modifications  related  to  the  new  input  require¬ 
ments  of  the  BLITS  model,  a  front-end  input  processor  was  needed  to  provide 
realistic  interpreted  data  to  the  BLITS  model.  This  input  engine  is  the  AIM 
processor.  In  addition  to  various  input  processing  routines,  AIM  was  originally 
designed  to  process  the  outputs  of  BLITS  into  comprehensible  and  compressed 
output  formats  that  would  follow  the  same  structure  as  previous  WAVES- 
developed  output  formats.  Thus,  other  WAVES  codes  could  be  adopted  that 
used  the  WAVES  FastVIEW  output  formats. 

The  FastVIEW  format  permitted  visualization  of  cloud  scenes  by  combining 
both  extinction  and  limiting  path  radiance  information  in  a  single  file.  Two  such 
visualization  algorithms  were  developed  that  were  compatible  with  outputs  from 
either  AIM  or  WAVES.  The  first  was  designed  as  an  orthographic  viewing 
technique  in  which  all  lines  of  sight  through  the  modeled  volume  were  treated 
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using  parallel  lines.  The  second  technique  produced  a  perspective  view  of  the 
modeled  volume,  where  a  mathematical  pinhole  camera  model  was  used  to 
construct  lines  of  sight,  all  converging  on  the  camera  position,  through  the 
pinhole  “aperture”  and  on  to  an  image  plane  within  the  system  (O’Brien  and 
Tofsted  1998).  To  construct  the  appropriate  input  data  sets  to  run  these  models, 
the  original  WAVES  output  data  sets  had  to  be  expanded  in  number.  Originally 
WAVES  produced  two  data  sets  comprising  a  gridded  surface  illumination  out¬ 
put  file  and  a  FastVIEW  output  file  of  gridded  3-D  extinction  and  directionally 
varying  limiting  path  radiance  information.  A  third  output  file  consisting  of  high- 
resolution  (62.5-m  resolution  vs.  the  standard  250-m  resolution)  extinction  data 
was  added  because,  when  visualized,  the  250-m  FastVIEW  data  appeared  as 
large  illuminated  blocks.  Using  the  higher-resolution  extinction  data,  along  with 
interpolation  methods  applied  to  the  coarser  resolution  FastVIEW  radiance  data, 
realistic  visualizations  of  cloud  features  and  illumination  could  be  produced. 

AIM  model  enhancements  for  AMIP 

The  goal  of  the  AMIP  program  from  the  perspective  of  the  AIM  models  was 
to  expand  the  AIM  modeling  capabilities  to  allow  surface  solar  loading  calcula¬ 
tions  as  a  user  option  and  to  permit  these  calculations  to  be  passed  along  to  other 
codes  in  the  form  of  solar  loading  databases.  These  databases  were  to  include 
information  on  direct  radiation  (unscattered  radiation  directly  from  the  solar  disk) 
as  well  as  diffusely  scattered  radiation.  The  goal  of  this  division  was  to  permit 
the  directionally  dependent  illumination  of  various  objects  within  a  terrain/ 
surface  model.  Using  this  information,  one  would  be  able  to  differentially 
illuminate  portions  of  the  terrain,  depending  on  surface  self-shadowing  features 
or  surface  directional  orientations  such  as  terrain  slopes  and/or  complex  vege¬ 
tation  surfaces.  The  data  sets  of  surface  illumination  created  are  2-D  in  nature, 
exhibiting  variations  in  position  due  to  overlying  clouds  and  their  shadows.  Thus, 
issues  such  as  the  time-dependent  changes  in  solar  loading  and  the  impact  on  the 
terrain  temperature,  and  how  these  variations  affect  target/background  signatures, 
can  be  addressed.  To  support  this  increasing  role  of  AIM,  a  series  of  changes 
and/or  additions  were  made  to  AIM.  The  following  sections  deal  with  the  nature 
and  implications  of  these  changes. 

Correlated  K-distribution  technique 

The  solar  radiation  that  provides  the  primary  heat  source  for  the  earth’s  sur¬ 
face  consists  of  visible  (0.35-0.75  pm)  and  near-infrared  (NIR)  (0.75 — 4.0  pm) 
energy.  Most  of  the  energy  at  the  shorter  wavelengths  (UV)  is  absorbed  by  ozone 
in  the  upper  atmosphere,  and  because  AIM  focuses  on  tropospheric  effects,  these 


Solar  Flux  Initialization  Schemes 


5 


absorption  effects  are  largely  accounted  for  by  initializing  the  AIM  models  with 
the  MODTRAN  code.  However,  across  the  NIR  spectrum  there  exist  several 
regions  where  water  vapor  absorption  is  important.  Concerns  over  propagation 
effects  across  these  band  structures  are  the  motivating  factor  in  moving  from  a 
series  of  sampled  absorption  effects  used  previously  in  AIM-visible  band  calcu¬ 
lations  to  the  integrated  solar  loading  calculations  needed  for  this  project. 

The  AIM  models  are  based  on  use  of  Beer’s  law,  which  assumes  monochro¬ 
matic  transmission  effects  (negative  exponential  function  of  the  optical  depth 
traversed).  Unfortunately,  when  considering  various  absorption  bands  in  the  NIR, 
the  optical  depth  can  vary  greatly  over  a  given  band.  Brute  force  methods  simply 
increase  the  number  of  spectral  channels  used  to  accommodate  the  increased 
variability  in  extinction  coefficients.  But  this  technique  is  possible  only  in 
models  that  make  concessions  in  accuracy  elsewhere,  such  as  using  a  lower 
spatial  dimensionality.  That  is,  the  model  atmospheres  are  not  allowed  to  vary 
horizontally.  To  permit  horizontal  variations  in  a  model  requires  commensurate 
increases  in  processing  time  because  the  computational  difficulty  is  proportional 
to  the  number  of  processing  cells  used.  The  250-m  resolution  typically  used 
within  BLITS/BLIRB,  though  rather  coarse  from  a  visualization  perspective, 
nevertheless  captures  typical  minimum  cloud  dimensions.  Yet  even  this  size  cell 
imposes  significant  computational  burdens. 

The  problem  was  how  to  permit  both  the  integration  over  the  full  solar 
spectrum  needed  to  compute  the  energy  arriving  at  the  surface  that  included 
horizontally  inhomogeneous  illumination  conditions  while  simultaneously 
avoiding  excessive  computational  times.  The  solution  to  these  competing 
priorities  came  in  the  form  of  an  adaptation  to  the  correlated  ^-distribution 
method  (Lads  and  Oinas  1991,  Fu  and  Liou  1992). 

The  correlated  ^-distribution  method  is  an  extrapolation  of  the  ^-distribution 
method.  Under  the  ^-distribution  (KD)  method,  instead  of  representing  the  varia¬ 
tions  in  the  extinction  coefficient  (k)  as  a  function  of  wavelength,  k  is  given  as  a 
function  of  its  probability  distribution,  the  probability  that  k  will  be  less  than  a 
given  value.  Under  these  circumstances,  the  governing  probability  parameter  g 
will  vary  between  0  and  1  as  Ovaries  between  minimum  and  maximum  values 
for  a  given  spectral  band  in  a  given  height  interval  above  the  surface.  The  major 
assumption  of  this  technique  is  then  invoked,  whereby  the  extinction  coefficient 
distribution  functions  at  different  levels  in  the  atmosphere  are  assumed  to  be  cor¬ 
related.  The  k  distributions  at  each  level  can  then  all  be  parameterized  through  a 
single  parameter  g.  Then,  instead  of  evaluating  the  propagation  characteristics  as 
a  function  of  wavelength,  the  same  quantities  can  be  evaluated  using  fewer  inter¬ 
vals  over  the  range  of  k  variations  as  functions  of  g. 
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We  took  this  correlation  technique  several  steps  further.  The  inputs  required 
to  evaluate  irradiance  values  within  the  AIM  program  BLITS  include  information 
concerning  the  diffuse  and  direct  radiation  at  the  upper  boundary  of  the  BLITS 
modeled  volume  as  well  as  layer  extinction  information.  Thus,  there  are  elements 
of  the  input  set  that  correspond  directly  to  data  produced  in  the  correlated  k- 
distribution  method  as  well  as  radiance  information.  But  we  would  expect  this 
radiance  data  to  be  somewhat  correlated  to  the  extinction  data  because  high- 
extinction  regions  should  correspond  with  low-radiance  regions  if  the  energy  is 
being  removed  in  overlying  layers  due  to  the  same  absorption  processes.  In  the 
typical  AJM  case,  the  modeled  volume  upper  boundary  is  at  4  km.  There  is  ample 
atmosphere  above  4  km  in  which  water  vapor  effects  absorb  sufficient  energy 
from  the  incident  radiation  to  produce  correlations  between  the  incident  radiation 
and  the  extinction  coefficient  data. 

To  assess  these  correlations,  we  had  to  alter  the  way  the  MODTRAN  model 
was  called  in  initializing  the  boundary  conditions  (at  the  4-km  boundary)  for  runs 
of  the  AIM  BLITS  model.  Originally,  for  the  visible  light  cases  studied,  MOD¬ 
TRAN  was  called  only  at  a  minimal  sampling  of  wavenumbers.  The  visible  band 
generally  spans  the  wavenumber  interval  between  14,440  and  24,240  cm-1.  The 
visible  band  was  modeled  as  seven  bands  of  1400  cm-1  by  sampling  the  spectral 
extinction  and  incident  radiance  at  the  ends  of  these  bands  using  a  100-cnf 1 
interval  and  averaging  the  maximum  values.  In  the  visible  band,  this  approach 
appeared  reasonable  due  to  the  low  absorptivity  of  the  atmosphere  over  each 
band.  But  in  the  calculation  of  the  solar  loading,  this  approach  would  be  inap¬ 
propriate  over  the  NIR  region  due  to  band  absorption  problems.  We  could  not 
greatly  increase  the  number  of  equivalent  bands  used  without  greatly  reducing 
performance,  nor  could  we  use  the  standard  correlated  ^-distribution  approach, 
which  primarily  applies  to  narrow  band  calculations.  In  fact,  because  our  run 
mechanism  relies  on  MODTRAN,  which  has  a  maximum  resolution  of  only  1 
cm-1,  it  is  impossible  to  model  the  actual  correlated  k  distribution  (CKD)  because 
that  calculation  requires  spectral  resolutions  on  the  order  of  anywhere  between 
20  and  2000  per  inverse  centimeter. 

Instead,  what  was  developed  was  an  approximation  to  the  correlated  k 
method.  In  this  approximation,  we  first  focused  on  transmission  rather  than  the 
extinction  coefficient  because  a  band  region  that  had  a  low  transmission  would 
be  similar  to  another  low-transmission  result  regardless  of  whether  the  extinction 
coefficients  were  off  by  an  order  of  magnitude  or  more.  The  significant  issue  was 
not  the  absolute  coefficient  value  but  rather  whether  energy  was  transmitted  or 
not,  since  we  were  primarily  interested  in  the  energy  reaching  the  surface. 

Second,  we  needed  to  consider  all  the  transmission  results  across  the  entire  solar 
band  at  once  because  we  had  to  reduce  the  number  of  channels  to  a  minimum. 
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Standard  CKD  methods  tended  to  use  narrow  bands  (perhaps  as  small  as  a  few 
hundred  wavenumbers  in  width),  however,  so  we  developed  measures  of  com¬ 
pactness  for  sets  of  sample  input  data  that  would  be  used  in  dividing  the  data 
between  the  different  channels. 

The  data  required  to  run  the  BLITS  model  include  an  equivalent  channel 
extinction  coefficient  for  all  molecular  absorbers  in  each  1-km  atmospheric  layer 
from  the  model  base  to  its  top.  In  addition,  values  for  the  incident  direct  irradi- 
ance  and  diffuse  radiances  in  a  set  of  downward-streaming  directions  are  needed 
for  boundary  input  information  at  the  volume  top.  In  general,  we  use  a  24-stream 
diffuse  radiance  model,  requiring  12  downward  diffuse  stream  calculations. 
When  combined  with  the  extinction  data  requirement  and  the  direct  radiance 
information,  14  MODTRAN  model  runs  in  all  were  required. 

Under  this  AMIP  program  we  looked  at  the  problem  of  model  initialization 
from  two  aspects.  The  first  considered  how  best  to  provide  data  needed  to  ini¬ 
tialize  the  thermal  model,  SWOETHERM,  in  order  to  provide  a  basis  for  model 
performance  analysis.  Considering  this  need,  we  wanted  to  determine  the  results 
for  several  case  studies  as  quickly  as  possible,  and  although  we  would  have 
preferred  to  develop  a  highly  optimized  code  for  rapid  computation,  it  was  not 
possible  given  the  timetable  required  to  provide  the  data.  Therefore,  we  initially 
developed  a  technique  for  making  the  CKD  computations  in  a  way  that  required 
as  few  model  alterations  as  possible.  As  a  result  of  performing  these  initial  cal¬ 
culations,  however,  we  discovered  that  the  processing  times  required  to  perform 
these  KD  computations  were  responsible  for  increasing  the  overall  model  com¬ 
putation  time  by  approximately  50%.  Thus,  after  the  case  study  data  had  been 
generated,  we  developed  a  second  method  for  evaluating  an  equivalent  set  of 
MODTRAN-based  boundary  and  initial  condition  data  that  will  permit  solar 
loading  calculations  but  which  would  not  require  the  direct  running  of  MOD¬ 
TRAN  at  execution  time.  We  have  divided  the  comments  regarding  code  im¬ 
provement  between  this  subsection,  which  only  discusses  the  original  CKD 
technique  itself,  and  additional  techniques  under  development  to  optimize  the 
processing  of  the  CKD  data. 

The  CKD  method  adapted  for  this  task  involved,  first,  the  development  of  a 
conceptual  model  of  the  input  data  required  to  run  the  radiative  transfer  model 
BLITS  over  the  entire  solar  spectrum.  In  that  model  it  was  first  realized  that 
correlating  14  independent  data  variables  (the  12  diffuse  stream  radiances,  the 
direct  stream  radiance,  and  the  four  1  -km  extinction  coefficients)  would  be  too 
much  to  attempt.  Further,  we  knew  that,  in  addition  to  the  data  produced  by  the 
MODTRAN  runs,  the  frequency  variable  was  also  of  interest  because  the  aerosol 
scattering  properties  in  the  troposphere  vary  strongly  with  wavelength.  We  thus 


8 


ERDC/CRREL  TR-03-13 


conceived  a  means  of  reducing  the  number  of  data  elements  to  be  correlated  from 
effectively  14  down  to  only  three.  The  method  chosen  was  to  interpret  a  single 
incident  energy  value  for  the  net  radiance  crossing  the  volume  upper  boundary  to 
represent  the  12  diffuse  channel  elements  plus  the  direct  channel  data.  It  was 
assumed  that  channel  data  undergoing  similar  extinction  or  scattering  would 
exhibit  similar  forms  for  the  net  directional  properties  of  the  scattered  energy. 

The  results  in  radiance  were  then  normalized  such  that  the  results  varied  between 
0  and  1 ,  representing  the  minimum  and  maximum  points  within  the  MODTRAN 
output  set.  A  second  variable  reflected  the  average  transmission  of  each  1-km 
layer.  A  third  variable  represented  a  normalized  wavenumber  parameter,  which 
varied  from  0  to  1  across  the  band  of  interest.  These  resulting  three-element  data 
sets  were  then  envisioned  as  pointing  to  locations  in  a  3-D  volume  where  wave- 
number  varied  along  one  axis,  the  generalized  incident  radiation  parameter  along 
a  second,  and  the  transmission  parameter  along  the  third. 

Running  the  MODTRAN  algorithm  at  2  cnf 1  resolution  over  the  solar  spec¬ 
trum  produced  a  “cloud”  of  points  within  a  unit  cube  once  the  results  for  each 
point  were  translated  into  their  associated  three  parameter  values  of  wavenumber, 
irradiance,  and  transmittance.  A  pattern  recognition  algorithm  was  then  used  to 
determine  logical  clusters  of  points  within  this  cloud  of  sample  values. 

The  method  of  running  the  pattern  recognizer  was  relatively  simple.  First,  all 
the  points  were  divided  into  a  minimal  number  of  “classes”  of  points  that  were 
generally  divided  into  wavebands  according  to  visible  color  classes  and  a  single 
NIR  band.  The  mean  values  of  each  of  the  three  parameters  were  then  computed, 
and  the  class  with  the  most  variance  about  the  means  was  divided  into  two 
daughter  classes  whose  mean  positions  differed  slightly  from  one  another  along 
the  axis  of  the  vector  component  with  the  greatest  variance.  Each  sample  point 
was  then  compared  to  the  newly  developed  class  centroids  and  new  class  assign¬ 
ments  were  determined,  depending  on  which  class  centroid  the  point  was  closest 
to.  Once  all  points  were  reassigned  to  a  new  class,  using  the  class  assignment 
value  of  each  point,  a  new  centroid  could  be  computed  for  each  class.  Subse¬ 
quently,  the  variance  of  points  within  each  class,  along  each  dimensional  axis, 
could  be  computed.  At  each  stage,  the  class  with  the  greatest  variance  in  any 
dimension  was  then  divided  into  two  new  classes,  and  points  belonging  to  the  old 
classes  would  be  reassigned  between  the  new  classes,  depending  on  their  proxim¬ 
ity  to  the  new  centroids.  Classes  would  continue  to  be  divided  until  one  of  two 
conditions  was  satisfied:  either  a  maximum  (user-specified)  number  of  classes 
had  been  reached,  or  (highly  unlikely)  all  variances  within  all  classes  reached 
zero.  In  the  case  of  our  computations,  since  we  would  never  be  able  to  process 
a  sufficient  number  of  classes  to  eliminate  all  variance  within  the  various  compu¬ 
tational  spectral  “channels,”  the  best  that  could  be  achieved  was  to  minimize  the 
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amount  of  variation  over  each  channel.  Overall,  this  pattern  recognition  tech¬ 
nique  is  a  variant  of  the  well-known  ISODATA  algorithm  (Ball  and  Hall  1 965). 

A  major  assumption  of  this  approach  was  that  variations  within  a  single 
channel  would  tend  to  cancel,  leading  to  overall  realistic  surface  solar  loading 
calculations,  even  though  results  for  individual  narrow-band  constituents  of  a 
channel  might  show  dispersion  about  that  result.  But  this  could  occur  only  if  the 
means  used  to  determine  the  average  statistics  for  the  channel  was  chosen  in  such 
a  way  as  to  be  connected  to  the  mean  associated  with  that  channel.  To  ensure  that 
class  statistics  were  situated  about  the  proper  mean  values,  the  wavelength  and 
transmittance  statistics  were  weighted  according  to  the  available  incident  energy 
at  each  individual  wavelength  interval.  For  example,  it  was  discovered  that  the 
transmission  and  incident  radiation  results  were  positively  correlated.  That  is, 
when  we  observed  energy  in  a  highly  absorbing  band  there  was  usually  less 
available  (more  had  already  been  absorbed  out  of  the  incident  energy  at  the  top 
of  the  modeled  volume)  than  in  a  nearby  low-absorption  region.  Thus,  it  would 
not  be  correct  to  simply  weight  the  transmission  parameter  evenly  across  a  single 
output  channel/class  of  data.  More  energy  would  be  available  for  propagation 
through  the  low-absorption  region,  and  thus,  overall  there  would  be  a  higher 
effective  transmission  capability  than  would  be  determined  if  properties  of  each 
wavelength  interval  were  given  equal  weight.  To  implement  this  technique 
within  AIM,  the  model  execution  stage  was  augmented  by  a  solar  loading  option, 
which  called  a  modified  version  of  MODTRAN  and  included  a  call  to  the  new 
classifier  algorithm. 

Y2K  compatibility  issue 

In  keeping  with  the  concern  over  Y2K  issues,  the  AIM  models  were  adjusted 
to  handle  crossover  into  the  next  century.  Previously,  these  codes  had  used  an 
ephemeris  computational  engine  that  was  destined  to  fail  for  dates  beyond  1999, 
but  this  failure  was  primarily  due  to  lunar  position  coefficients  that  were  based 
on  a  1977  effective  date.  Subsequently  these  codes  were  updated  to  permit 
unlimited  computations  into  the  future  based  on  recomputation  of  the  governing 
coefficients  for  the  given  year  of  the  run  being  made.  There  remain,  however, 
possible  differences  between  the  algorithms  used  by  AIM  and  those  used  under 
SWOETHERM.  The  differences  have  not  been  completely  accounted  for,  but 
appear  small  relative  to  solar  loading  concerns.  They  may  be  related  to  observer 
altitude  issues. 
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SWOETHERM  model  interoperability 

To  build  in  compatibility  between  the  AIM  outputs  and  the  SWOETHERM 
inputs,  we  added  an  auxiliary  routine  to  translate  between  the  standard  250-m 
resolution  surface  information  output  file  produced  by  AIM  and  a  standard  1 00- 
m  format  expected  by  SWOE.  To  provide  this  compatibility,  we  had  the  option 
of  increasing  the  resolution  of  the  AIM/BLITS  radiative  transfer  model  run  or 
reformatting  the  output.  We  chose  to  reformat  the  data  rather  than  run  at  higher 
resolution.  This  was  a  practical  rather  than  an  optimal  decision.  Using  a  100-m 
resolution  grid  within  BLITS  would  have  produced  output  databases  1 5.6  times 
larger  than  the  current  data  sets,  which  typically  occupy  300  MB  of  disk  space. 
Such  an  expansion  seemed  unrealistic. 

Instead,  a  cosine  interpolation  routine  previously  used  for  generating  image 
overlays  was  adapted  to  produce  100-m-resolution  solar  loading  maps.  These 
files  contained  information  about  direct,  diffuse,  and  surface-emitted  total 
radiation  over  the  wavelength  interval  of  interest.  Of  course,  since  there  is  little 
thermal  radiation  over  the  solar  band,  the  emitted  values  were  uniformly  near 
zero.  In  addition,  these  files  contained  information  about  the  horizontal  dimen¬ 
sions  of  each  cell  and  the  overall  domain  of  the  modeled  volume  (standardized 
at  8  km  in  both  the  east-west  and  north-south  directions).  In  addition,  the  refor¬ 
matting  technique  transformed  the  order  in  which  the  results  were  output.  The 
standard  AIM  output  presented  the  data  in  columns  running  south  to  north  as  the 
inner  loop.  The  SWOETHERM  format  required  the  data  to  be  presented  in  terms 
of  rows  of  data  running  east  to  west  as  the  inner  loop  and  with  the  columns  pre¬ 
sented  in  a  north-to-south  order. 

AIM  model  acceleration 

As  mentioned  above,  in  extending  the  AIM  model  to  evaluate  effects  over 
the  entire  solar  spectrum,  MODTRAN  had  to  be  run  at  rather  high  spectral 
resolution  (2  cnT1)  and  the  results  post-processed  using  the  pattern  recognition 
classifier  based  on  correlated  ^-distribution  theory.  The  unfortunate  side  effect 
of  this  procedure  was  that  while  we  were  able  to  provide  AIM  with  the  inputs  it 
needed  to  run  the  BLITS  model,  the  “cure”  was  rather  costly.  Running  MOD¬ 
TRAN  at  high  resolution  took  nearly  as  long  as  the  combined  running  times  of 
the  remainder  of  the  AIM  codes.  Thus,  from  a  practical  standpoint,  the  continued 
use  of  MODTRAN  becomes  problematic  to  the  use  of  these  modules  in  any 
realistic  real-time  or  other  Army  operational  application.  However,  if  that  pur¬ 
pose  is  ultimately  to  be  achieved,  it  will  be  necessary  to  disconnect  AIM  or 
any  future  3-D  radiative  processing  techniques  from  MODTRAN  or  other  high- 
spectral-resolution  processing  routine.  The  issue  is  this:  high-spectral-resolution 
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codes  assume  a  1-D  atmosphere  that  is  simple  to  model  and  uses  a  coarse  spatial 
resolution.  The  product  of  the  spatial  elements  times  the  number  of  spectral 
channels  is  then  roughly  comparable  for  AIM  and  MODTRAN,  but  while  we 
have  control  of  AIM,  permitting  us  to  take  shortcuts  on  the  spatial  side,  MOD¬ 
TRAN  must  be  preprocessed  to  eliminate  its  processing  costs  up  front.  While 
these  costs  are  then  translated  into  the  storage  space  required  to  archive  the 
resulting  data  sets,  at  least  they  would  not  impact  processing  time  for  the  actual 
computations. 

The  overall  goal  of  this  computational  engine  is  to  generate  any  output  pro¬ 
duct  within  a  maximum  of  5  minutes.  Although  we  are  not  going  to  get  to  this 
point  at  the  conclusion  of  this  program,  we  are  considerably  closer  than  we  were 
previously  and,  with  the  advent  of  faster  computers  with  more  memory,  the  pros¬ 
pect  of  achieving  such  performance  or  better  within  the  next  10  years  is  highly 
likely.  What  follows  then  are  the  steps  that  were  achieved  during  this  program  in 
terms  of  model  speed  increases. 

The  first  performance  enhancements  were  applied  to  die  BLITS  algorithm. 
Though  these  changes  tended  to  slightly  decrease  the  model  accuracy  (1-3% 
greater  error),  the  resulting  model  produces  results  that  are  still  well  below  the 
uncertainties  in  several  other  model  input  variables  such  as  cloud  cover,  cloud 
properties,  haze  layer  thickness,  and  haze  extinction  properties.  Two  changes 
involved  aspects  of  the  BLITS  main  routine.  First,  the  error-checking  coefficients 
were  increased  significantly.  Previously  a  point  value  difference  between  two 
iterations  of  0.01%  was  considered  sufficient  to  flag  an  item  as  an  error.  This 
condition  criterion  was  increased  by  a  factor  of  20,  reducing  the  number  of  model 
iterations  required  to  satisfy  stability  of  the  solution  by  approximately  30%. 

The  second  major  change  focused  on  a  means  of  accelerating  the  algorithm 
convergence.  Previously,  BLITS  simply  iterated  on  the  scattering  source  using 
the  previous  iteration’s  diffuse  radiance  terms  as  input  to  the  current  iteration’s 
scattering  equations  and  continued  this  iteration  process  until  the  diffuse  results 
stabilized.  However,  this  method  was  improved  upon  by  using  an  acceleration 
technique  operating  on  the  results  of  two  consecutive  iterations.  This  technique 
requires  additional  memory  to  store  a  new  copy  of  the  diffuse  stream  array  from 
an  intermediate  computation,  but  it  permits  the  iteration  procedure  to  converge 
with  a  third  less  iteration. 

A  final  BLITS  acceleration  technique  regarded  a  change  in  the  means  of 
initializing  the  effects  of  direct  beam  illumination  on  individual  cells  within  the 
modeled  volume.  Previously,  lines  of  sight  were  traced  to  sample  points  along 
the  outer  walls  of  each  cell  within  the  modeled  volume.  Statistical  averages  were 
then  generated  to  describe  overall  illumination  for  each  cell  and  for  the  overall 
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direct  energy  scattered,  absorbed,  and  transmitted  within  each  cell.  However,  in 
cases  where  the  sun  was  near  the  horizon,  this  technique  became  computationally 
intensive  due  to  the  horizontal  periodicity  normally  employed  in  describing  the 
edge  effects.  Under  these  conditions,  for  each  point  on  an  input  cell  surface  a  line 
of  sight  to  the  volume  upper  boundary  could  span  many  copies  of  the  modeled 
volume  prior  to  exiting  the  volume. 

To  accelerate  this  process,  the  results  at  each  height  level  in  the  volume 
would  be  stored  and  used  at  the  next  lower  level  in  the  ray  tracing  process.  For 
each  new  cell,  interpolations  on  results  for  previously  evaluated  cells  could  be 
used  so  that  a  line  of  sight  need  go  through  no  more  than  several  cells  at  most.  In 
many  cases,  results  could  be  obtained  by  passing  only  through  a  single  cell.  This 
change  led  to  slightly  longer  times  for  cases  in  which  the  sun  was  directly  over¬ 
head,  but  much  faster  processing  times  for  cases  in  which  the  sun  was  near  the 
horizon.  Given  that  the  sun  must  be  near  the  horizon  at  least  twice  each  day,  this 
change  saves  considerable  processing  time  on  average.  There  is  a  slight  increase 
in  the  model  error  due  to  the  cascading  effects  of  interpolations  on  interpolations 
at  each  lower  level.  In  some  cases,  these  cascading  effects  are  minimized  (for 
small  zenith  angles),  but  for  other  cases  [large  zenith  angles  with  low  density  (no 
cloud)  atmospheres]  the  errors  tend  to  be  on  the  order  of  several  percent. 

Eliminating  direct  MODTRAN  runs 

The  final  acceleration  method  is  dictated  by  the  need  to  avoid  direct  MOD¬ 
TRAN  runs  in  calculating  solar  loading.  The  reason  for  this  need  is  obvious: 
MODTRAN  must  be  run  at  high  resolution  to  obtain  information  on  the  NIR 
absorption  band  structure  that  is  critical  in  determining  the  total  energy  reaching 
the  surface.  However,  this  same  accuracy  requirement  is  a  liability  to  rapid  cal¬ 
culation  of  effects.  Thus,  direct  calls  to  MODTRAN  should  be  avoided.  To 
eliminate  them,  it  will  be  necessary  to  produce  sets  of  preprocessed  data 
representing  characteristic  MODTRAN  outputs  under  a  variety  of  conditions. 
Then  interpolations  can  be  made  on  the  precomputed  data  sets. 

The  first  step  in  this  process  was  to  understand  the  nature  of  the  data  space  of 
MODTRAN  outputs.  Then  a  method  of  interpolation  was  developed  that  was  ex¬ 
tensible,  given  a  variable-size  output  database.  First,  consider  the  size  of  the  data¬ 
base.  Tests  were  conducted  where  MODTRAN  was  exercised  under  a  varying  set 
of  conditions.  Output  variables  of  interest  were  those  required  by  AIM  as  input: 
extinction  coefficients,  direct  irradiance,  and  diffuse  radiance  at  the  top  of  the 
AIM  modeled  volume.  For  purposes  of  these  computations,  AIM’s  volume  was 
fixed  at  a  4-km  vertical  extent.  Initially  only  solar  source  data  sets  were  consid¬ 
ered.  Presumably,  if  the  model’s  upper  boundary  was  allowed  to  vary  in  height 
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and  a  lunar  source  was  permitted,  then  this  would  extend  the  database  require¬ 
ment  by  two  additional  dimensions.  As  it  was,  the  extinction  coefficients  were 
found  to  vary  only  with  respect  to  the  vertical  structure  model  selected  (1976 
U.S.  standard  atmosphere,  midlatitude  summer,  etc.)  and  height  of  the  surface 
above  sea  level.  The  solar  direct  radiation  was  found  to  vary  with  solar  zenith 
angle  in  addition  to  the  original  two  variables.  The  diffuse  parameters  showed 
additional  dependence  on  the  surface  reflectivity  and  haze  layer  visibility. 

Four  of  these  data  variables  could  be  evaluated  as  floating-point  numbers  for 
interpolation  purposes  (zenith  angle,  surface  reflectivity,  visibility,  and  surface 
height  above  sea  level).  But  the  atmosphere-type  parameter  was  provided  to 
MODTRAN  in  terms  of  an  integer.  However,  it  was  noticed  that  these  “model” 
values  associated  with  various  atmospheric  state  scenarios  (midlatitude  summer, 
midlatitude  winter,  subarctic  summer,  subarctic  winter,  U.S.  standard  atmos¬ 
phere,  and  tropical)  could  be  described  according  to  latitude  and  season,  which 
could  be  expressed  in  terms  of  two  additional  parameters.  The  six  main  atmos¬ 
pheric  scenarios  could  be  associated  with  2-D  season/latitude  data  points  and 
interpolations  made  with  respect  to  a  user-specified  point  within  the  grid.  In  one 
sense,  then,  these  results  would  permit  a  user  greater  flexibility  than  was  current¬ 
ly  possible  using  the  existing  MODTRAN  code,  because  what-if  questions  could 
be  asked  about  intermediate  conditions,  while  the  previous  model  permitted  only 
fixed  choices  of  atmospheric  models  to  be  used. 

It  was  decided  that  the  resolution  in  each  of  the  data  variables  was  a  signifi¬ 
cant  consideration  for  interpolation  purposes.  That  is,  having  too  few  cases  could 
result  in  bad  interpolations.  Also,  some  variables  had  more  dynamic  range  than 
others.  For  example,  zenith  angle  could  vary  between  0  and  90°,  but  surface 
reflectivity  could  vary  only  between  0  and  1 .  To  account  for  these  dynamics, 
different  variables  were  resolved  by  different  numbers  of  sample  points.  Zenith 
was  varied  between  0  and  88°  in  increments  of  1 1°.  Visibility  was  varied  from  1 
km  through  64  km  in  logarithmic  fashion  (2  km,  4  km,  8  km,  etc.).  Reflectivity 
was  varied  linearly  between  0.05  and  0.20  in  increments  of  0.05;  it  was  then  set 
to  0.3, 0.5,  and  0.9  to  handle  cases  of  highly  reflective  surfaces  (white  sand 
beaches  and/or  snow-covered  surfaces)  albeit  at  lower  resolution.  The  six  unique 
model  cases  enumerated  above  were  each  computed  separately,  but  categorized 
according  to  a  scheme  that  assigns  tropical  results  to  a  latitude  of  10°,  midlatitude 
results  to  latitude  40°,  and  subarctic  results  to  latitude  70°.  Since  the  tropical 
results  are  presumed  not  to  be  affected  by  seasonal  variation,  they  are  assigned 
to  the  tropical  latitude  regardless  of  season.  The  summer  results,  however,  are 
assigned  to  the  Julian  date  associated  with  the  summer  solstice,  and  similarly  the 
wintertime  results  are  assigned  to  the  winter  solstice.  The  U.S.  standard  atmos¬ 
phere  is  then  treated  as  a  special  case  and  assigned  to  the  midlatitude  equinox 
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condition.  Finally,  the  height  above  sea  level  of  the  model  surface  was  ranged 
from  0  to  2  km  in  increments  of  250  m. 

Overall,  there  are  nine  zenith  positions,  seven  visibilities,  seven  reflectivities, 
six  models,  and  nine  altitudes,  for  a  total  of  23,814  cases  to  be  computed  under 
this  archival  approach.  Obviously,  it  is  not  practical  to  archive  the  complete 
MODTRAN  outputs  of  all  these  runs,  nor  would  it  be  wise  to  await  the  comple¬ 
tion  of  computing  results  for  every  case  before  beginning  to  attempt  interpola¬ 
tions  on  a  partial  data  set.  It  is  estimated  that  the  total  processing  time  required 
to  evaluate  this  23,814-case  data  set  is  over  one  year  on  a  typical  workstation. 
Instead,  first,  the  CKD-processed  channel  data  for  a  given  case  had  to  be  com¬ 
puted  and  the  resulting  data  stored  for  this  case.  And  second,  a  method  of  inter¬ 
polation  had  to  be  developed  that  could  produce  results  for  partial  databases. 

Let  us  consider  the  first  aspect  of  this  problem.  Because  we  are  interpolating 
our  results  between  difference  conditions,  for  which  CKD  results  must  be 
archived,  the  CKD  method  described  previously,  which  optimizes  the  choice 
of  channel  assignment  for  each  2  cm-1  interval,  cannot  be  used  as-is  in  this 
archiving  mode.  The  reason  is  the  channel  assignment  technique.  It  would  be 
highly  likely  (and  it  was  noted  in  tests  of  the  method)  that  similar  (though  not 
exact)  conditions  produced  variable  channel  assignments.  Thus  any  interpolation 
method  used  must  rely  on  a  means  of  archiving  the  MODTRAN  results  that 
regularizes  the  method  of  channel  assignment.  To  do  this,  a  separate  version  of 
the  classifier  algorithm  was  developed  that  read  the  choice  of  channel  assign¬ 
ments  from  an  external  file.  This  file  was  from  a  run  of  the  original  model  on 
what  was  considered  a  standard  case:  average  solar  zenith  angle,  average  surface 
reflectivity,  etc.  The  assignment  data  for  this  case  were  output  and  stored  in  a 
separate  file  that  was  then  used  as  the  template  for  channel  assignments  in  all 
subsequent  classifier  calculations. 

Because  it  might  be  possible  to  run  the  surface  radiation  calculations  at 
varying  resolution,  die  results  of  the  MODTRAN  runs  were  processed  using 
variable  numbers  of  channels:  First,  results  were  computed  for  the  original  seven- 
channel  visible  scenario,  then  solar  loading  calculations  were  made  with  7, 15, 
20,  and  25  channels.  Each  of  these  five  scenarios  is  created  using  the  output  from 
MODTRAN  runs  over  the  solar  band  for  every  combination  of  the  main  five 
variables  described  above. 

Interpolating  results  using  sample  MODTRAN  runs 

The  major  challenge  associated  with  this  technique  is  to  perform  meaningful 
interpolations  based  on  less  than  a  full  set  of  sample  case  outputs.  The  problem  is 
simple:  given  a  total  of  23,814  cases  and  a  production  rate  of  approximately  65 
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cases  per  24-hr  period  under  optimal  conditions,  the  database  will  require  just 
over  a  year  of  continual  processing  time  to  evaluate,  assuming  no  stoppages  in 
production.  What  if  the  computational  engine  is  halted  temporarily?  How  can  the 
order  of  computation  be  optimized  such  that  interpolations  can  be  as  accurate  as 
possible  as  soon  as  possible? 

To  accommodate  computational  interruptions,  a  simple  technique  was  devel¬ 
oped  whereby  a  file  listing  all  previously  computed  results  is  updated  each  time 
a  new  set  is  added.  The  interpolation  method  also  refers  to  this  listing  in  order  to 
interpolate  input  conditions  using  results  computed  to  date. 

To  optimize  the  order  of  the  computations,  a  two-stage  process  was  used.  In 
the  first  stage,  48  data  sets  were  calculated,  which  represented  the  outer  shell  of 
the  most  extreme  conditions  to  be  modeled.  For  the  zenith  angle,  the  0  and  88° 
cases  were  considered.  This  meant  that,  for  the  first  four  data  elements,  minimum 
and  maximum  conditions  were  considered.  For  the  atmosphere-type  parameter, 
the  subarctic  winter  and  summer  conditions  were  considered,  along  with  the 
tropical  condition,  for  a  total  of  three  conditions,  which  mapped  into  five  data 
sets  (the  tropical  case  maps  into  three  cases  total  due  to  seasonal  independence). 
The  48  computed  cases  thus  map  into  a  total  of  80  cases. 

Once  this  initial  series  of  extreme  cases  were  run,  the  processing  code  con¬ 
sidered  interior  points  using  a  weighting  scheme.  For  each  index  position  of  each 
variable,  a  number  between  0  and  4  was  assigned,  indicating  a  somewhat  arbi¬ 
trary  assessment  of  the  importance  of  evaluating  results  at  that  value  of  the 
variable  in  question.  For  example,  for  the  ground  altitude  parameter,  the  0  m 
above  sea  level  condition  was  given  highest  priority  (weight  4),  followed  by  the 
1000-m  condition  (weight  3),  and  then  by  the  500-m  condition  (weight  2),  etc. 
The  weights  for  each  combination  of  run  conditions  are  then  added  and  a  sum  is 
obtained.  The  maximum  sum  of  all  the  weights  is  20.  The  order  in  which  compu¬ 
tations  are  performed  is  determined  by  an  iterative  procedure.  First,  only  condi¬ 
tions  whose  weight  total  sums  to  20  are  computed.  Once  all  the  possible  combi¬ 
nations  have  been  checked  and  all  the  20  weight  results  have  been  computed,  the 
algorithm  cycles  back  again  and  computes  all  the  cases  whose  weight  total  sums 
to  19.  Following  this  iterative  approach,  presumably,  the  most  important  cases 
will  be  computed  first.  In  this  way  the  algorithm  gradually  fills  in  the  set  of 
computed  results  in  a  smooth  manner  and  avoids  the  problems  associated  with 
systematically  processing  the  model  space  in  a  standard  iterative  fashion  that 
might  leave  large  data  gaps  until  late  in  the  production  process. 

Given  a  database  in  a  state  of  partial  completion,  what  method  of  interpola¬ 
tion  is  best?  To  answer  this  question,  we  initially  considered  full  interpolation 
over  an  array  of  nearby  points,  but  found  that  this  technique  would  not  be  able  to 
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account  for  gaps  in  the  available  data.  The  problem  is  balance:  given  a  partially 
filled  database,  it  is  nearly  impossible  to  build  a  window  correctly  around  a  point 
of  interest  that  equally  balances  the  influence  of  sample  points  on  opposite  sides 
of  the  point  in  question.  Equally  important  are  methods  of  dealing  with  points 
near  the  boundaries  of  the  hypercube  of  possibilities. 

For  such  cases  it  is  important  to  develop  a  method  that  will  not  merely 
average  results  from  nearby  points.  To  perform  an  interpolation,  let  us  assume 
the  vector  x  -  is  composed  of  vector  components  xy,  where  i  cycles  from  1  to  6, 
representing  the  six  different  parameters  (zenith  angle,  visibility,  etc.),  and  where 
j  cycles  from  1  to  J,  the  number  of  sample  data  sets  currently  processed.  As  time 
increases,  J  will  increase,  but  at  processing  time  this  number  will  be  fixed. 

Associated  with  each  point  in  the  parameter  space  there  would  be  a  specific 
vector  output  of  MODTRAN  that  describes  the  input  needed  to  run  AIM  at  a  set 
of  spectral  channels.  Let  the  elements  of  this  vector  (Y )  be  given  as  components 
(Yk).  Thus  we  have  a  space  where  components  are  functions  of  the  point  in  the 
parameter  space,  or. 


?=/(*),  c 

The  problem  is  then  to  solve  for  Y  given  results  Y  ■  at  a  finite  set  of  sample 
points  Xj .  The  question  is  how  best  to  accomplish  this  task. 

Originally  an  approach  was  considered  whereby  a  quadratic  interpolation 
about  a  given  point  could  be  computed.  However,  the  solution  of  this  problem 
would  involve  high-order  terms  that  could  not  be  readily  solved  for,  and  the 
result  could  be  unstable.  Instead,  a  model  was  developed  that  was  linear  in  each 
of  the  six  parameters  and  centered  about  the  point  in  question.  That  is,  we 
assume, 


^  (*)  _  Ar,o  +  X \,i (x;  Xj).  (2) 

/=o 

X  is  the  point  in  parameter  space  of  interest,  with  the  Xt  being  the  com¬ 
ponents  of  X  .  The  A  coefficients  are  produced  as  a  result  of  a  least  squares 
analysis.  There  will  be  a  separate  set  (.4  j  of  coefficients  for  each  element  YK 
of  the  solution  vector. 

By  centering  the  solution  about  point  X  ,  once  the  A  coefficients  have 
been  solved  for  through  the  least  squares  analysis,  and  since  we  are  interested 
in  the  point  x  =  X ,  we  then  have  the  simple  solution,  YK  (x)  =  AK  fi .  Thus  the 
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formality  of  the  least  squares  process  is  simply  to  remove  any  linear  trends  in  the 
data  such  that  the  resulting  estimate  of  the  value  will  not  contain  any  remaining 
linear  bias. 

We  solve  for  the  Arf  using  the  available  data  Y  (xj )  J .  However,  we  also 
need  to  account  for  the  nature  of  the  data  set  available.  A  normal  least  squares 
analysis  would  give  each  sample  point  equal  weight.  And  while  we  would  as¬ 
sume  that  the  functional  dependence  of  Y  should  be  relatively  slowly  varying 
with  any  of  the  individual  parameters,  it  may  be  that  the  overall  dependence  with 
some  of  the  parameters  may  be  nonlinear.  In  this  case,  giving  equal  weight  to  all 
sample  points  would  be  bad,  because  distant  points  would  tend  to  reduce  the 
estimate  (for  negative  second  derivative)  or  raise  the  estimate  (for  positive 
second  derivative).  Thus,  it  is  necessary  to  weight  the  value  of  available  data 
points  such  that  linearities  in  the  nearby  region  would  receive  greater  weight. 

For  example,  let  the  distance  between  the  central  point  (z)  and  a  sample 
point  (*.)  be  measured  as 


i  now  varies  between  0  (the  constant  term)  and  6.  Now  let  the  weighting  scheme 
assume  the  form 


where  Wj  is  the  weight  assigned  to  the  y'th  sample  and  a  is  a  width  parameter 
that  is  fitted  to  the  distribution  of  computed  data  sets  about  the  point  of  interest, 

X  .  Given  the  weighting  scheme  above,  let  us  define  the  sum-squared  error  as 

Ar 1  (5) 

j= 1  L  1=0 

where  AK  is  the  error  in  the  Afh  element  of  Y . 

From  a  theoretical  point  of  view,  this  approach  assumes  that  each  parameter 
is  independent  of  all  others.  Obviously  this  leads  to  calculations  that  make  no 
sense,  such  as  solar  zenith  angles  of  0°  at  subarctic  latitudes,  but  as  long  as  cr  is 
kept  sufficiently  small,  this  should  not  lead  to  significant  errors. 

Using  the  least  squares  analysis  approach,  we  take  the  partial  derivative  of 
Ak  with  respect  to  one  of  the  variables  Ay , .  The  resulting  problem  requires  the 
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solution  of  a  matrix  multiplication,  B  =  MA  ,  where  M  is  a  square  matrix  com¬ 
posed  of  correlation  elements. 

_ Although  the  solution  of  this  matrix  problem,  involving  the  matrix  inverse  of 

M ,  is  fairly  routine,  there  remains  the  issue  of  determining  the  width  parameter, 
<j ,  given  a  particular  arrangement  of  nearby  sample  points.  If  a  is  set  to  too 
large  a  value,  faraway  points  may  have  too  significant  an  influence  on  the  final 
results.  For  example,  if  there  are  numerous  nearby  points  already  available,  it 
makes  little  sense  to  permit  distant  points  to  have  much  bearing  on  the  outcome. 
On  the  other  hand,  if  there  are  few  nearby  points,  it  may  be  necessary  to  increase 
a  in  order  to  ensure  a  statistical  sample  that  is  meaningful.  The  question  is,  how 
should  cr  be  determined? 

Compounding  this  problem  is  the  fact  that  we  are  not  dealing  with  a  “nor¬ 
mal”  three-dimensional  space.  Rather,  the  interpolation  space  is  six-dimensional 
in  nature.  And  since  a  is  used  to  weight  the  total  distance  in  any  combination  of 
dimensions,  we  need  to  know  the  hypervolume  of  a  hypersphere  of  radius  a  in 
6-D  space  to  be  able  to  determine  how  dense  the  set  of  sample  points  is  in  the 
neighborhood  of  X  .  To  determine  the  functional  form  of  this  space,  a  point¬ 
counting  algorithm  was  developed.  Using  this  method  it  was  determined  that  a 
sphere  of  radius  r  in  six  dimensions  has  a  hypervolume  of  V6(r)  =  ttV6  /  6 .  This 
result  was  used  in  an  algorithm  that  measures  the  density  of  available  points  in 
the  vicinity  of  the  point  of  interest  [X  J  for  hyperspheres  of  different  radii. 

If  there  is  a  full  set  of  available  points  within  a  certain  range,  then  the  sample 
density  will  be  near  unity.  If,  however,  there  is  a  scarcity  of  nearby  points,  the 
density  at  small  r  values  will  be  much  less  than  1 .  It  will  then  be  necessary  to 
determine  how  the  density  increases  to  a  maximum  at  some  larger  range.  As  long 
as  the  set  of  sample  data  is  less  than  complete,  this  density  must  fall  below  unity 
at  some  range  of  r  values.  It  therefore  makes  sense  to  choose  a  as  the  minimum 
distance  at  which  the  ratio  F(r)  =  N(r)/V6(r)  is  maximized,  where  N(r)  is  the 
number  of  sample  points  within  range  r.  That  is,  we  choose  cr  to  be  the  mini¬ 
mum  range  r  at  which  the  density,  F(r) ,  is  a  local  maximum.  Of  course,  if  we  set 
r  large  enough,  V6{r)  will  always  be  much  larger  than  N(r)  because  the  array  of 
computed  results  will  always  be  finite  in  nature,  whereas  V6(r)  can  increase  to 
infinity.  Thus  we  expect  that  there  will  always  be  some  local  maximum,  but  that 
F{r)  may  remain  relatively  close  to  its  maximum  value  for  some  range  of  dis¬ 
tances  as  the  array  becomes  more  fully  populated  with  computed  results.  It 
therefore  makes  sense  to  choose  the  minimum  value  at  which  this  maximum  is 
reached  so  that  points  at  large  distances  are  not  overweighted  in  the  least-squares 
method. 
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SWOE  field  program  calculations 

To  demonstrate  model  interoperability  between  the  AIM  models  and  the 
SWOE  thermal  model,  meteorological  information  collected  during  the  three 
SWOE  field  programs  (Grayling  I,  Yuma,  and  Grayling  II)  were  used  in 
generating  AIM  surface  solar  loading  data  sets  for  selected  canonical  weather 
conditions.  Three  cloud  scenarios  (clear,  partly  cloudy,  and  overcast)  were 
chosen  for  each  data  site.  Hourly  surface  weather  observations  were  used  to 
determine  primary  model  inputs.  In  addition,  a  morning  rawinsonde  data  set  was 
used  in  estimating  cloud  layer  thickness  and  base  height.  We  calculated  the  direct 
and  diffuse  components  of  the  solar  loading  for  an  8-  by  8-km  grid  at  a  0.25-km 
spatial  resolution  and  then  used  a  smoothing  interpolation  method  to  produce  a 
0.1 -km  resolution  output  grid  compatible  with  the  SWOE  models. 

Semi-empirical  surface  radiation  calculations 

Shapiro  (1972, 1982, 1987)  has  developed  a  simple  model  to  determine  the 
direct  and  diffuse  shortwave  fluxes  at  the  surface  using  only  standard  surface 
meteorological  observations.  This  model  has  been  used  successfully  to  initialize 
surface  energy  budget  models.  Unlike  the  AIM  model,  this  model  does  not  pro¬ 
duce  spatially  distributed  solar  fluxes.  For  clear  skies  and  overcast  conditions, 
using  nondistributed  flux  values  may  have  little  impact  on  distributed  surface 
energy  budget  models,  provided  the  spatial  domain  is  not  too  large.  Using  non¬ 
distributed  fluxes  for  partly  cloudy  conditions  raises  some  interesting  issues. 

Model  description 

The  basic  model  approach  involves  dividing  the  atmosphere  into  k  layers  and 
assuming 

Rk  +  Tk  +  Ak  =  \  (6) 

where  R,  T,  and  A  are  the  reflectance,  transmission,  and  absorption  of  layer  k, 
respectively.  R,  T,  and  A  have  been  parameterized  in  terms  of  the  solar  zenith 
angle  60  and  the  state  of  the  atmosphere/clouds  using  the  very  extensive  SOL- 
MET  database  (NOAA  SOLMET  Vol.  2  1979).  The  general  form  of  the  flux 
equations  is 


'i=vir1+vi  <7a> 


(7b) 
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Layer  k  =  0  is  the  top  of  the  atmosphere.  For  operational  use,  the  atmosphere  has 
been  divided  into  three  layers  consistent  with  the  concept  of  low,  middle,  and 
high  clouds,  where  layer  1  corresponds  to  the  high  cloud  region,  layer  2  to  the 
middle  cloud  region,  layer  3  to  the  low  cloud  region,  and  layer  4  is  the  ground. 
The  downwelling  solar  flux  at  the  ground  (bottom  of  layer  3)  is  given  as 

Ia=T,T2T,/D2  (8) 

and 

D2  =d3(dtd2  -RiR2T^)-dlR2apidT^  - R^jTjJ  (9) 

and 


dj  - 1  RjRj+] 


(10) 


a ^  is  the  ground  albedo.  Assuming  the  total  transmission  Tk  can  be  specified  as 
the  sum  of  the  transmission  of  the  direct  solar  component  7tir  and  the  diffuse 
solar  component  T^. 


Tk  =  Tkdir +  Tff 


The  direct  solar  flux  component  at  the  ground  is  given  as 


(11) 


j  dir  _  rpdir  >jidir  rj^dir  j 

12  3 


and  the  diffuse  component  as 


If 


(12) 


(13) 


Ioi  (=1369.3  W/m2)  is  the  shortwave  flux  at  the  top  of  the  atmosphere.  To  solve 
the  above  equations  for  the  downwelling  direct  and  diffuse  flux,  it  is  necessary  to 
assume 


T*f  =  Rk  (14) 

and  therefore  Tkr  can  be  obtained  from 

rj^dir  _  rj-i  _ J^dif 


(15) 
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Shapiro  (1972, 1982, 1987)  parameterized  Tk  and  Rk  in  terms  of  the  atmos¬ 
pheric  and  cloud  conditions  as  follows: 


& 

II 

£ 

+ 

►—* 

1 

iL 

** 

(16) 

Tk  =  <PkTk+(l-<pk)tk 

(17) 

i? 

II 

-V 

(18) 

pk  is  the  cloud  reflectance  for  cloud  layer  k,  rk  the  clear  sky  reflectance  for  layer 
k,  vk  the  cloud  transmission  for  layer  k,  tk  the  clear  sky  transmission  for  layer  k, 
fk  the  fractional  cloud  amount  for  layer  k,  and  IT  is  a  cloud  weighting  factor,  p, 
r,  t,  r ,  and  W  are  parameterized  in  terms  of  the  cosine  of  the  solar  zenith  angle 
using  the  SOLMET  data  set. 

Pk  =  ak  +  al  cos  0o  +  ak  cos2  +  al cos3  9o 

(19a) 

rk=a°k  +  a\  cos  90  +  a\  cos2  90  +  a\  cos3  9a 

(19b) 

tk=a°k  +  a\  cos 90+a2k  cos2  90  +  a\  cos3  90 

(19c) 

?k  =  al  +  a\  cos  9o  +  al  cos2  9a  +  a\  cos3  90 . 

(19d) 

ak',  aak,  bk,  and  bbk  are  parameterized  in  terms  of  the  following  atmospheric  and 
cloud  categories:  clear,  smoke  and  haze,  thin  cirrus  and  cirrostratus,  thick  cirrus 
and  cirrostratus,  altostratus  and  altocumulus,  and  low  clouds.  The  cloud  weight¬ 
ing  factor  W  is  given  as 

W  =  c0+clcos90  +  c2fk  +  cjk cos 9  0  +  c4 cos2  6 0  +  cjk  .  (20) 

The  cs  are  parameterized  in  terms  of  the  following  cloud  categories:  thin  cirrus 
and  cirrostratus,  thick  cirrus  and  cirrostratus,  altostratus  and  altocumulus,  and 
low  clouds.  The  value  of  the  coefficients  can  be  found  in  Shapiro  (1987). 
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Atmospheric  transmission  model:  MODTRAN 

MODTRAN  is  a  MODerate-resolution  TRANsmittance  (Berk  et  al.  1989) 
model  for  predicting  atmospheric  radiance  and  transmittance.  Specifically,  it 
calculates  atmospheric  transmittance,  single-scattered  solar  and  lunar  radiance, 
direct  solar  and  lunar  irradiance,  and  multiple-scattered  solar  and  lunar  radiance. 
MODTRAN  calculates  continuum-type  and  line  (utilizing  band  models  that  are  a 
function  of  pressure,  temperature,  and  line  width)  molecular  absorption,  molecu¬ 
lar  scattering,  and  aerosol  and  hydrometer  absorption  and  scattering.  The  spectral 
resolution  of  the  model  is  2  cm-1  (full  width  at  half-maximum)  in  averaged  steps 
of  1  cm-1.  The  spectral  fluxes  calculated  by  MODTRAN  are  integrated  over  the 
appropriate  spectral  interval  to  generate  the  broadband  solar  fluxes.  MODTRAN 
is  basically  a  plane  parallel  model  and  cannot  compute  the  fluxes  for  partly 
cloudy  skies  directly.  This  is  achieved  by  weighting  the  MODTRAN  fluxes 
for  clear  and  overcast  sky  conditions  with  the  appropriate  cloud  amounts. 

Measured  solar  flux  values 

The  measured  fluxes  have  been  obtained  from  the  SWOE  data  set.  Under  the 
SWOE  program,  field  campaigns  were  conducted  at  two  locations  for  three  sepa¬ 
rate  time  periods.  The  Grayling  site  chosen  to  represent  the  NATO  European 
analog  is  on  one  of  the  tank  firing  ranges  in  the  Camp  Grayling  Military  Reserva¬ 
tion,  northeast  of  Grayling,  Michigan.  This  area  is  typically  rolling  hills  with  a 
mixture  of  vegetation  types  from  bare  soils  to  grasses  to  forests.  The  underlying 
soils  are  predominantly  deep  sands  with  little  organic  material  in  the  near-surface 
layers.  The  specific  area  used  in  the  SWOE  Grayling  I  (September  and  October 
1992)  and  Grayling  II  (4  March  1994  to  15  April  1994)  field  programs  is  a  valley 
oriented  roughly  northeast-southwest,  characterized  by  an  open,  somewhat  grass- 
covered  floor.  Ridges  to  the  east  and  west  are  predominantly  covered  by  mixed- 
growth  trees  (mostly  pine  and  oak).  This  valley  is  about  2  km  long  by  1  km  wide. 
Except  for  very  low  sun  elevation  angles,  the  surrounding  ridges  do  not  shadow 
the  pyrheliometer  and  pyranometer  used  to  collect  the  solar  flux  information. 

The  site  at  Yuma  Proving  Ground,  Arizona,  chosen  for  SWOE’s  southwest 
Asian  analog,  was  characterized  as  upland  Sonoran  desert  with  primary  and 
secondary  washes  throughout.  Soils  were  varied  from  highly  graded  coarse 
gravel  in  the  wash  areas  to  large  expanses  of  desert  pavement  on  undisturbed 
ridges.  Underlying  soils  were  predominantly  mixed  gravel  and  sands  over  the  test 
area.  Vegetation  ranged  from  short  grasses  to  cacti  to  large  palo  verde  trees.  The 
wash  pattern  was  evident  due  to  the  abundance  of  vegetation  in  these  low  areas. 
At  the  start  of  the  test,  the  entire  region  was  extremely  wet  after  several  months 
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of  higher  than  average  rainfall.  During  the  later  part  of  the  field  program,  it  was 
very  dry  and  during  the  day  there  was  considerable  airborne  dust. 
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3  THERMAL  (SWOETHERM)  MODEL 

SWOETHERM  is  a  one-dimensional  state-of-the-ground  model  that  is  an 
aggregate  of  three  interactive  models.  SWOETHERM  consists  of  a  soil/snow 
thermal  model,  a  vegetation  (grasses)  thermal  model,  and  a  canopy  (deciduous  or 
coniferous)  thermal  model.  When  vegetation  or  canopy  exists  the  modifications 
of  fluxes  due  to  the  interaction  of  the  vegetation  or  canopy  with  the  underlying 
soil  surface  are  taken  into  consideration. 

Soil/snow  thermal  model 

The  soil/snow  model  used  in  SWOETHERM  is  SNTHERM  (Jordan  1991), 
developed  by  Rachel  Jordan.  SNTHERM  is  a  one-dimensional  mass  and  energy 
balance  model  for  predicting  the  temperature  profile  within  strata  of  snow  and 
soil.  It  considers  the  transport  of  liquid  water  and  water  vapor  and  the  phase 
changes  of  water  as  components  of  the  heat  balance  equation.  The  impacts  of 
snow  accumulation,  ablation,  densification,  and  metamorphosis  on  the  snow 
thermal  and  optical  properties  are  modeled.  The  infiltration  of  water  in  the  snow 
is  modeled  assuming  gravity  flow.  When  snow  is  present,  the  water  infiltrating  to 
the  snow/soil  interface  is  artificially  drained  from  the  system.  The  snow  and  soil 
are  divided  into  horizontally  infinite  control  volumes,  and  the  mass  and  heat 
balance  equations  are  applied  to  each  control  volume.  A  spatial  discretization 
scheme  similar  to  a  finite-difference  method  is  used  in  the  spatial  domain,  while 
a  Crank-Nicolson  method  is  used  to  discretize  the  time  domain.  The  model  uses 
an  adaptive  time-step  procedure  that  automatically  adjusts  the  time  step  (typically 
between  900  and  5  seconds)  to  obtain  the  desired  accuracy  of  the  solution  to  the 
mass  and  heat  balance  equations.  The  governing  equations  are  linearized,  and  a 
tridiagonal-matrix  algorithm  is  used  to  obtain  the  desired  solution.  The  required 
meteorological  boundary  conditions,  which  are  either  user-supplied  or  generated 
internally  in  the  model,  are  air  temperature,  dew  point  temperature,  wind  speed, 
precipitation,  and  the  solar  and  downwelling  infrared  fluxes.  The  soil  and  snow 
layer  physical  properties  required  by  the  model  can  either  be  supplied  by  the  user 
or  will  default  to  the  hardwired  values  in  the  code.  The  same  is  true  for  the  pro¬ 
files  of  temperature  and  water  in  the  snow  and  soil  layers. 

Vegetation  thermal  model 

The  energy  budget  of  a  simple  vegetation  layer  on  a  soil  surface  is  modeled 
using  a  steady-state  semi-infinite  plane  parallel  model  (Balick  et  al.  1981),  which 
is  described  by  the  foliage  emissivity  and  albedo,  a  foliage  height,  the  foliage 
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fractional  coverage  and  a  foliage  state  parameter.  The  vegetation  consists  of  a 
single  homogeneous  layer  that  is  infinite  in  both  x  and  y  directions.  The  sum  of 
the  energy  terms,  consisting  of  the  absorbed  solar  and  infrared  fluxes,  the  emitted 
longwave  flux,  and  the  sensible  and  latent  heat  fluxes,  is  equal  to  0  at  each  time 
increment.  The  solution  of  the  resulting  polynomial  equation  of  degree  n  for  the 
foliage  temperature  is  obtained  using  a  root-finding  algorithm  or  by  linearization 
of  the  equation.  The  original  model  used  a  root-finding  algorithm,  but  more 
recently  it  has  been  modified  to  use  a  linearization  procedure.  The  foliage  frac¬ 
tional  coverage  and  the  shortwave  albedo  are  a  function  of  the  vegetation  type 
(high,  medium,  or  low)  and  the  season  (winter,  spring,  summer,  or  fall)  and  are 
hardwired  in  the  software. 

Canopy  thermal  model 

The  canopy  model  (Smith  et  al.  1981)  is  a  semi-infinite,  steady-state  plane 
parallel  energy  budget  model  consisting  of  three  canopy  layers,  an  atmospheric 
layer  above  the  canopy,  and  a  ground  layer  below  the  canopy.  The  model  con¬ 
siders  the  longwave  and  shortwave  fluxes  and  the  interactions  between  the 
various  layers  and  uses  a  simplified  expression  for  the  sensible  heat  and  evapo- 
transpiration  flux  calculations.  The  model  also  incorporates  the  orientation  and 
distribution  of  leaves  in  the  canopy  layers.  Both  a  longwave  and  a  shortwave 
transfer  matrix  are  computed  for  each  canopy  layer  based  on  a  leaf  frequency 
distribution  model  and  the  leaf  area  index.  This  matrix  is  used  to  compute  the 
transfer  of  longwave  and  shortwave  fluxes  within  the  canopy  and  the  interaction 
of  the  canopy-layer-emitted  longwave  flux. 
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4  ANALYSIS 


Instrument  variability 

To  compare  the  spatial  variations  of  the  measured  and  AIM-derived  solar 
flux  values  for  Grayling  and  Yuma,  it  is  necessary  to  remove,  or  at  least  quantify, 
variations  due  to  inconsistencies  between  the  instruments.  The  ten  Eppley  8-48s 
and  five  Eppley  PSPs  used  to  obtain  the  solar  flux  measurements  during  the 
Yuma  test  were  located  at  a  single  site,  and  data  were  collected  under  clear-sky 
conditions.  Figure  1  is  a  plot  of  the  maximum,  mean,  minimum,  and  standard 
deviation  for  each  instrument  type  for  the  Yuma  measurements.  At  approximate¬ 
ly  0800  and  1730  local,  the  standard  deviation  and  the  measured  solar  flux  values 
show  a  relatively  large  spike.  This  is  an  unknown  artifact  and  may  be  associated 
with  sun  glint  off  part  of  the  instrument.  The  standard  deviation  for  both  instru¬ 
ment  types  is  on  the  order  of  5  to  10  W/m2,  but  there  is  approximately  a  70-W/m2 
difference  between  the  mean  of  the  two  instrument  types  around  solar  noon.  The 
maximum  difference  between  the  minimum  and  maximum  flux  values  for  both 
instrument  types  was  approximately  28  W/m2  and  occurred  around  1230  local. 
The  percent  difference  between  either  the  maximum  or  minimum  and  the  mean 
between  0900  and  1700  local  was  less  than  3.5%. 

A  similar  analysis  was  performed  for  instrumentation  used  during  the  Gray¬ 
ling  II  experiment.  Only  a  single  instrument  type  was  used  at  the  Grayling  II  field 
program.  The  measurements  were  made  using  14  PSPs  taking  measurements  over 
an  approximately  8-hr  period.  The  data  prior  to  0900  local  vary  considerably  (see 
Figure  2).  Some  of  this  variation  may  be  due  to  shadowing  associated  with  a  low 
ridgeline  to  the  east  of  the  site.  After  approximately  1030  local,  the  standard 
deviation  remains  fairly  constant,  showing  only  a  slight  increase  with  time  from 
about  10  W/m2  to  about  13  W/m2.  The  maximum  difference  between  the  maxi¬ 
mum  and  minimum  flux  values  was  approximately  55  W/m2  and  occurred  around 
local  noon.  The  percent  difference  between  the  maximum  and  the  mean  or  the 
minimum  and  the  mean  after  approximately  0930  local  was  on  the  order  of  5%. 
Data  measured  before  0930  local  is  unreliable. 
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For  Yuma,  the  standard  deviation  and  the  maximum  difference  for  each  set 
of  instrument  types  are  similar.  Somewhat  unexpected  is  the  rather  large  mean 
difference  between  the  two  instrument  types.  This  problem  will  not  be  encoun¬ 
tered  for  Grayling  because  only  a  single  instrument  type  was  used.  The  maxi¬ 
mum  difference  for  the  Grayling  data  set  is  about  twice  the  maximum  difference 
for  the  Yuma  data  set.  In  general,  it  can  be  anticipated  that  5%,  and  possibly  as 
high  as  10%,  of  the  variance  in  the  solar  flux  measurements  at  any  specific  time 
can  be  attributed  to  instrument  differences.  One  word  of  caution:  the  percent 
differences  and  standard  deviations  are  fairly  large  for  solar  zenith  angles  greater 
than  80°. 

Solar  flux  analysis 

To  assess  the  utility  of  the  three  solar  flux  calculation  methods,  canonical 
cloud  data  sets  were  selected  from  each  of  the  three  SWOE  field  programs.  Each 
model  produces  direct,  diffuse,  and  total  solar  flux  values,  but  it  is  primarily  the 
total  solar  fluxes  that  are  compared.  Shapiro’s  technique  generates  fluxes  based 
on  a  semi-empirical  technique.  The  resulting  fluxes  depend  on  the  cloud  condi¬ 
tions,  visibility,  geographical  location,  and  the  time  of  year  and  day.  MODTRAN 
is  a  plane-parallel  model  that  depends  on  the  model  atmosphere,  visibility  and 
aerosol  type,  location,  and  time  of  year  and  day.  Partly  cloudy  conditions  are 
taken  into  consideration  by  running  the  model  for  clear  and  overcast  conditions 
and  weighting  the  results  according  to  the  actual  cloud  amount.  AIM  produces 
spatially  distributed  solar  fluxes  for  each  time  step,  while  the  other  models  and 
the  measured  values  consist  of  a  single  flux  value  for  each  time  step.  The  SWOE- 
THERM  model  was  run  for  bare  ground,  grass,  snow,  and  deciduous  and  conif¬ 
erous  canopies  using  either  the  measured  solar  fluxes  or  the  model-calculated 
solar  fluxes.  Only  the  cloud  cover  varied  from  model  run  to  model  run.  All  other 
meteorological  conditions  were  held  constant.  The  model  was  run  for  a  24-hr 
period  starting  at  midnight  for  each  surface  condition.  The  resulting  temperature 
values  for  each  surface  condition  and  cloud  cover  are  compared. 

Clear  sky  flux  analysis 

For  the  clear-sky  cases,  the  AIM  solar  flux  values  were  constant  over  the 
extent  of  the  spatial  domain  (8  km  by  8  km).  The  other  models  and  the  measured 
data  also  consist  of  a  single  solar  flux  value. 
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Grayling  I  -  25  Sept  92  (DOY  269) 
Total  Solar  -  Clear 


Figure  3.  Clear-sky  flux,  both  measured  and  modeled,  based  on  meteoro¬ 
logical  conditions  from  Grayling  I  field  program. 


Grayling  I:  Clear  skies 

Figure  3  represents  the  total  solar  flux  for  clear-sky  conditions  for  Grayling  I. 
It  is  believed  that  the  dip  in  the  measured  solar  flux  values  between  1400  and 
1600  local  is  associated  with  either  thin  cirrus  or  a  haze  layer.  The  observer  re¬ 
ported  a  haze  layer  aloft  after  1 700  local.  The  maximum  difference  between  the 
measured  total  solar  flux  and  the  total  flux  calculated  using  Shapiro’s  technique 
is  only  13  W/m2,  and  the  average  difference  is  only  -2.3  W/m2. 

A  comparison  for  the  same  time  period  of  two  meteorological  sites  separated 
by  several  hundred  meters  indicates  a  maximum  difference  of  25  W/m2  and  an 
average  difference  of  10  W/m2.  Both  the  AIM-calculated  total  flux  and  the 
MODTRAN  total  flux  values  are  considerably  less  than  the  measured  values. 

The  maximum  difference  and  average  differences  between  the  measured  and  the 
AIM-calculated  fluxes  are  1 14  W/m2  and  43  W/m2,  respectively.  For  the  MOD¬ 
TRAN  fluxes,  the  maximum  and  average  differences  are  1 1 1  W/m2  and  26  W/m2, 
respectively.  Although  it  is  not  entirely  clear  why  there  should  be  differences  of 
this  magnitude,  the  most  likely  explanation  focuses  on  how  the  models  handle 
the  aerosol  loading  in  the  lower  levels  of  the  atmosphere.  The  Shapiro  model 
does  not  characterize  the  surface  haze  visibility  at  all,  while  both  the  MODTRAN 
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and  AIM  models  account  for  varying  haze  layer  optical  depth  via  standardized 
vertical  profiles  based  on  input  surface  layer  visibility.  To  explore  the  relation¬ 
ship  between  input  visibility  and  solar  loading,  we  ran  the  AIM  model  using  its 
standard  2-km-thick  haze  layer  under  varying  visibility  conditions.  In  Figure  4, 
as  the  visibility  increases  (smaller  optical  depths),  the  direct  and  total  solar  com¬ 
ponents  increase,  while  the  diffuse  component  increases  until  the  visibility  ex¬ 
ceeds  3  km.  The  decrease  in  the  diffuse  component  for  visibilities  greater  than  3 
km  is  due  to  a  decrease  in  the  scattering  of  the  direct  component.  The  solar  flux 
reaching  the  surface  is  sensitive  to  the  visibility,  especially  when  the  visibility  is 
less  than  20  km.  During  most  of  the  day  selected  for  this  comparison,  the  average 
visibility  was  around  28  km. 


Figure  4.  Change  in  solar  fiux  components  as  function  of  visibility  as  com¬ 
puted  using  AIM  model. 

Table  1  summarizes  the  surface  temperature  difference  between  the  measured 
and  modeled  solar  fluxes  for  different  surface  materials.  The  measured  minus 
Shapiro  maximum  temperature  difference  (Maxdif)  does  not  exceed  0.5°  for  any 
material  type.  The  average  difference  (Avgdif)  is  even  less.  The  relative  accuracy 
of  the  thermal  model  is  on  the  order  of  1°C.  The  measured  minus  MODTRAN 
and  measured  minus  AIM  maximum  temperature  difference  occurs  for  bare 
ground  and  is  on  the  order  of  1  to  2°  C  and  the  average  difference  is  on  the  order 
of  a  half  of  a  degree.  When  the  surface  temperature  differences  displayed  in 
Table  1  are  positive,  it  indicates  the  temperature  calculated  using  the  model  flux 
is  cooler  than  the  value  calculated  using  the  measured  flux.  For  example,  the 
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surface  temperature  differences  associated  with  the  measured  minus  the  AIM 
fluxes  are  all  positive,  indicating  the  AIM  mean  flux  values  are  always  less  than 
the  measured  flux  values.  For  this  clear-sky  case,  the  Shapiro  calculated  flux 
values  are  the  closest  to  the  measured  values  for  Grayling  I.  The  Shapiro  model 
is  also  the  model  that  has  the  smallest  computational  burden.  Thus,  for  clear  skies 
for  fairly  small  spatial  domains  and  for  relatively  high  visibilities,  this  model 
offers  a  rapid  method  of  computing  the  direct  and  diffuse  solar  flux  for  initial¬ 
izing  thermal  models. 


Table  1.  Surface  temperature  differences  (°C)  for  different  surface  types  and 
different  solar  flux  initialization  schemes  for  clear-sky  conditions  for  Grayling  1. 

Surface  conditions 

Measured 

Bare 

Snow 

Deciduous 

Coniferous 

Grass 

1.0 

0.7 

msM 

0.8 

-0.4 

-0.3 

-0.4 

MODTRAN 

0.5 

0.3 

0.3 

■91 

0.2 

0.1 

0.2 

mm 

0.2 

0.1 

-0.5 

-0.3 

-0.1 

-0.2 

Shapiro 

-0.1 

-0.1 

0.1 

0.1 

0.0 

mm 

mm 

1.0 

'mom 

1.0 

0.0 

WEm 

AIM 

i.i 

0.7 

0.5 

mm 

wEM 

Yuma:  Clear  skies 

Both  the  MODTRAN  and  Shapiro  total  solar  flux  values  are  greater  than  the 
measured  flux  values  and  the  AIM  total  solar  flux  values  are  less  than  the  corre¬ 
sponding  measured  value.  The  maximum  difference  between  the  measured  minus 
the  Shapiro  derived  flux  is  -45  W/m2,  and  the  average  over  the  daylight  period  is 
-20.5  W/m2.  The  maximum  difference  between  the  measured  minus  the  MOD¬ 
TRAN  flux  values  is  55.8  W/m2  and  occurs  at  0700  local.  All  other  MODTRAN 
values  for  the  time  period  are  greater  than  the  measured  values,  with  an  average 
difference  of -1 1 .2  W/m2.  The  maximum  and  average  difference  for  the  meas¬ 
ured  minus  the  AIM  total  solar  flux  values  is  57.4  and  30.4  W/m2,  respectively. 
The  maximum  measured  difference  for  two  sites  separated  by  several  hundred 
meters  was  28  W/m2;  the  average  difference  between  the  two  sites  was  4.0  W/m2. 
Again,  the  difference  between  measured  and  modeled  fluxes  is  greater  than  be¬ 
tween  two  measurement  sites.  In  general,  the  visibility  for  the  Yuma  clear  sky 
day  is  approximately  twice  that  for  the  Grayling  I  clear  sky  condition.  The 
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impact  of  the  different  solar  flux  values  for  different  surface  types  is  given  in 
Table  2.  As  anticipated,  the  maximum  temperature  difference  occurs  for  bare 
ground  and  is  on  the  order  of  the  relative  accuracy  of  the  thermal  model. 


Table  2.  Surface  temperature  differences  (°C)  for  different  surface  types 
and  different  solar  flux  initialization  schemes  for  clear-sky  conditions  for 
Yuma. 

Surface  conditions 

Measured 

Bare 

Grass 

1.4 

0.6 

■ 

-1.1 

-0.3 

MODTRAN 

-0.1 

0.0 

■ 

0.0 

■ 

-0.7 

Shapiro 

| 

-0.1 

Maxdif 

1.6 

1.2 

Mindif 

0.0 

-0.1 

AIM 

Avgdif 

0.4 

0.2 

Grayling  II:  Clear  skies 

The  maximum  differences  in  the  total  solar  flux  for  the  measured  minus 
MODTRAN,  measured  minus  Shapiro,  and  measured  minus  AIM  are  48.0  W/m2, 
15.5  W/m2,  and  85.4  W/m2,  while  the  average  differences  are  8.4  W/m2, 1 .8 
W/m2,  and  26.0  W/m2,  respectively.  The  average  and  maximum  difference  in  the 
total  solar  flux  for  two  meteorological  sites  separated  by  approximately  1 00  m 
for  the  same  time  period  was  10.0  W/m2  and  101.0  W/m2.  The  maximum  differ¬ 
ence  occurred  between  0725  and  0805  local  and  may  be  due  to  shadowing  of  one 
of  the  sites  by  a  low  ridgeline  to  the  east.  A  plot  of  the  calculated  SWOETHERM 
surface  temperatures  for  bare  soil  using  the  different  flux  initialization  schemes  is 
given  in  Figure  5. 

Again,  the  maximum  temperature  difference  and  the  greatest  average  differ¬ 
ence  occur  for  bare  soil  (Table  3).  In  general,  the  differences  are  not  greater  than 
the  relative  accuracy  of  the  thermal  model.  For  all  clear-sky  scenarios,  the  flux 
model  used  to  initialize  the  thermal  model  does  not  significantly  change  the  re¬ 
sulting  surface  temperature.  For  this  relatively  small  spatial  domain,  AIM  com¬ 
putes  a  single  value  for  the  solar  flux  for  each  time  step.  This  would  not  be  the 
case  for  relatively  large  domains  where  the  solar  elevation  angle  would  change 
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with  location.  For  clear  skies  and  relatively  small  spatial  domains,  a  simple  flux 
model  like  Shapiro’s  has  the  advantage  of  being  relatively  fast  and  accurate.  For 
larger  spatial  domains,  the  issue  of  changes  in  the  solar  zenith  angle  must  be 
addressed.  Even  in  these  cases,  a  model  like  Shapiro’s  could  be  run  for  several 
locations  over  the  domain  to  model  the  effects  of  solar  zenith  angle  variations  on 
the  solar  flux  at  the  surface. 


Grayling  li  -  Bare  ground,  clear  skies  - 16  Mar  94  (DOY  075) 


Hour  of  Day 

Figure  5.  Calculated  surface  temperature  values  for  bare  ground  and  differ¬ 
ent  solar  initialization  schemes. 


Cloudy  sky  flux  analysis 

To  investigate  the  impact  of  overcast  sky  conditions  on  the  calculated  solar 
flux  and  surface  temperatures,  days  with  low  overcast  cloud  conditions  were 
selected  when  possible.  The  direct  solar  flux  was  generally  less  than  100  W/m2 
and  frequently  0.  For  cloudy  conditions,  the  AIM  fluxes  varied  over  the  spatial 
extent  of  the  domain.  AIM  used  a  250-m  grid  resolution  that  was  interpolated  to 
100-m  resolution,  resulting  in  6400  flux  values  for  the  8-  by  8-km  domain  for 
each  1-hr  time  step.  Because  it  would  be  computationally  prohibitive  to  run 
SWOETHERM  for  all  6400  flux  values  for  each  hour,  the  following  procedure 
was  used.  For  each  hour,  the  maximum  and  minimum  flux  values  were  found, 
and  the  average  of  the  6400  values  was  determined.  Three  SWOETHERM  model 
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runs  were  then  made  for  each  hour:  one  using  the  maximum  value,  one  using  the 
minimum  value,  and  one  using  the  average  value. 


Table  3.  Surface  temperature  differences  (°C)  for  different  surface  types  and  dif¬ 
ferent  solar  flux  initialization  schemes  for  clear-sky  conditions  for  Grayling  II. 

Surface  conditions 

Measured 

Bare 

Snow 

Deciduous 

Coniferous 

Grass 

Maxdif 

0.3 

0.3 

0.2 

0.1 

0.2 

Mindif 

-0.4 

-0.3 

-0.1 

-0.1 

-0.2 

MODTRAN 

Avgdif 

0.1 

0.1 

0.0 

0.0 

0.0 

Maxdif 

0.2 

0.1 

0.1 

0.1 

0.1 

Mindif 

—0.3 

-0.2 

-0.1 

0.0 

-0.2 

Shapiro 

Avgdif 

0.0 

0.0 

0.0 

0.0 

0.0 

Maxdif 

0.9 

0.5 

0.5 

0.3 

0.5 

Mindif 

0.0 

0.0 

-0.1 

0.0 

0.0 

AIM 

Avgdif 

0.3 

0.2 

0.1 

0.0 

0.2 

Grayling  I:  Cloudy  skies 

The  distribution  of  the  AIM  total  solar  flux  for  1300  local  is  depicted  in 
Figure  6.  Solar  flux  measurements  were  made  at  several  sites  during  the  Grayling 
I  field  program.  All  sites  are  within  a  circle  with  a  radius  of  approximately  1  km. 
Thus  it  is  not  possible  to  determine  the  spatial  distribution  of  solar  flux  values 
from  these  sites  for  the  8-  by  8-km  area.  Instead,  we  have  equated  the  temporal 
variations  of  the  measured  fluxes  to  the  spatial  variations  as  computed  with  the 
AIM  model.  The  AIM  model  in  essence  represents  an  instant  in  time.  To  generate 
the  temporal  variations,  we  used  the  1-min  flux  values  from  a  single  site  for  a 
period  from  30  minutes  before  to  30  minutes  after  the  hour.  Based  on  the  winds 
at  the  cloud  height  and  assuming  only  cloud  advection  and  no  dynamics,  during 
this  1-hr  period  a  cloud  would  advect  approximately  20  km.  Using  this  concept, 
it  would  take  only  22  minutes  for  a  cloud  to  be  advected  across  the  region.  It  was 
felt  that  using  only  22  one-minute  observations  would  not  provide  a  sufficient 
sample  for  comparison  of  the  observed  data  and  the  AIM-calculated  values. 
Figure  7  is  a  histogram  of  the  6400  AIM  values  produced  for  conditions  at  1300 
local,  compared  to  the  distribution  of  the  60  observed  values.  The  AIM  model 
has  values  in  the  500  to  800  W/m2  range.  At  1300  local,  the  cloud  conditions 
were  0.9  stratocumulus  at  300  m  and  0.1  altocumulus  at  7000  m.  While  the  AIM- 
computed  flux  values  appear  higher  than  the  observed  conditions,  there  are 
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limitations  to  the  available  data  set.  First,  the  cloud  spatial  distribution  was  gen¬ 
erated  using  CSSM.  It  appears  that  CSSM  generated  areas  where  there  were  no 
clouds  or  the  clouds  were  created  with  small  optical  depths.  One  explanation  for 
lower  optical  depths  would  be  the  lack  of  information  on  the  actual  thickness  of 
the  cloud  layers.  Estimates  of  the  layer  thickness  were  developed  through  assess¬ 
ment  of  the  morning  rawinsonde  flight  data,  but  direct  measurements  were 
unavailable.  Another  explanation  is  the  way  cloud  fraction  information  is 
reported.  In  the  data  sets  provided,  the  total  fractional  cloud  cover  of  all  three 
(low,  medium,  and  high)  cloud  layers  never  exceeded  unity.  Thus,  higher  cloud 
layers  may  have  had  larger  fractional  coverage,  but  were  not  reported.  This 
would  result  in  an  erroneous  cloud  thickness.  Any  of  these  influences  would 
account  for  the  higher  flux  values,  and  points  to  a  need  for  improved  reporting 
requirements  during  future  testing.  Yet  even  for  the  measured  flux  values  at  a 
single  site,  the  fluxes  varied  by  a  factor  of  two  over  the  60-min  period,  indicating 
that  even  under  overcast  conditions  the  solar  flux  can  vary.  These  variations  may 
be  due  to  variations  in  the  cloud  optical  depth. 


Total  Solar  Radiation 

Grayling  1-11  Oct  92  (DOY  285)  1300  hr  -  Cloudy  - 
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Figure  6.  Distribution  of  total  solar  flux  at  1300  local  as  calculated  with  AIM 
model. 
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Figure  7.  Histogram  of  AIM  fluxes  for  1300  local  and  measured  fluxes  (see 
text  for  details). 

The  total  solar  flux  values  for  the  different  solar  radiation  initialization 
schemes  are  presented  in  Figure  8.  As  indicated  earlier,  at  approximately  1300 
local  the  reported  sky  conditions  consisted  of  nine-tenths  low  cloud  and  one- 
tenth  middle  cloud.  Figure  9  presents  the  resulting  surface  temperatures  from 
the  different  initialization  schemes  for  bare  soil.  The  AIM  initialization  scheme 
includes  the  use  of  the  mean,  maximum,  and  minimum  flux  values  for  each  time 
period.  In  general,  the  maximum  difference  in  the  surface  temperature  for  the 
MODTRAN,  Shapiro,  and  AIM  average  value  initialization  scheme  is  less  than 
2°.  This  is  not  the  case  for  the  initializations  using  the  AIM  maximum  and  mini¬ 
mum  fluxes.  The  maximum  values  result  in  a  surface  temperature  difference  that 
is  on  the  order  of  4°  warmer  while  the  minimum  is  approximately  4°  cooler  than 
the  temperatures  computed  using  the  measured  flux  values.  At  1300  local,  the 
difference  in  the  calculated  bare  soil  surface  temperature  using  the  maximum  and 
minimum  flux  is  approximately  8°.  This  represents  the  extreme  since  the  mini¬ 
mum  and  maximum  flux  values  are  used  for  the  entire  period.  That  is,  since  the 
maximum  and  minimum  are  used  for  the  entire  period  leading  up  to  1 300  local 
there  will  be  an  accumulative  effect.  In  reality,  a  single  location  in  the  simulation 
area  would  experience  a  range  of  flux  values  as  cloud  conditions  over  the  loca¬ 
tion  change.  We  would  expect  the  surface  temperature  to  vary  between  the  upper 
and  lower  limits  given  by  the  upper  and  lower  curves  in  Figure  9.  The  impact  of 


Solar  Flux  Initialization  Schemes 


37 


the  different  initialization  schemes  is  summarized  in  Table  4.  Negative  values 
imply  the  surface  temperature  computed  using  the  measured  flux  values.  Maxdif, 
Mindif,  and  Avgdif  are  the  maximum  and  minimum  temperature  difference  and 
the  average  of  the  differences  for  the  period  when  solar  fluxes  are  not  zero  (day¬ 
light).  As  anticipated,  the  greatest  differences  occur  when  we  initialize  the  model 
using  either  the  AIM  maximum  or  minimum  flux  values.  In  addition,  the  greatest 
differences  occur  for  bare  surface  conditions. 

Grayling  I  -  1 1  Oct  92  (DOY  285) 

Total  Solar  -  Cloudy 


Figure  8.  Total  solar  fluxes  for  overcast  conditions. 

With  regard  to  these  computations,  several  items  should  be  noted.  First,  the 
AIM  average  calculations  show  an  average  difference  that  is  consistently  closer 
to  the  measured  conditions  (within  approximately  0.1°)  than  either  the  MOD- 
TRAN  or  Shapiro  methods.  The  AIM  maximum  and  minimum  cases,  on  the 
other  hand,  tend  to  bracket  the  observed  conditions  equally  on  the  high  and  low 
side.  By  comparison,  the  Shapiro  predictions  consistently  produce  predicted  tem¬ 
peratures  above  those  predicted  using  the  measured  fluxes.  The  advantage  of  the 
Shapiro  model  is  that  it  requires  fewer  input  parameters  and  is  computationally 
fast.  The  MODTRAN  estimates  are  consistently  low,  probably  due  to  the  inabil¬ 
ity  of  MODTRAN  to  model  clouds  with  varying  optical  depths,  as  compared  to 
both  the  Shapiro  and  AIM  methods. 


38 


ERDC/CRREL  TR-03-13 


Grayling  -  Bare  ground,  cloudy  skies  - 11  Oct  92  (DOY  285) 


Figure  9.  Surface  temperature  for  bare  soil  calculated  using  different  initial¬ 
ization  schemes. 


Yuma:  Cloudy  skies 

During  the  day  the  cloud  conditions  varied  from  a  solid  low  overcast  with 
rain  and  drizzle  to  scattered  low  clouds  with  overcast  to  broken  mid-level  cloud 
cover.  Comparing  the  flux  values  (Figure  10),  it  appears  that  the  calculated  AIM 
average  flux  value  is  not  consistent  with  the  peak  associated  with  the  measured 
value  at  1000  local.  This  is  somewhat  misleading,  however,  since  the  variable 
plotted  is  the  average  flux  of  all  6400  values  over  the  8-  by  8-km  region.  In 
general,  we  are  comparing  point  measurements  and  point  model-generated  fluxes 
(MODTRAN  and  Shapiro)  with  spatially  distributed  model-generated  flux  values 
at  a  single  instant  of  time.  Figure  1 1  is  a  plot  of  the  maximum  flux  values  for 
each  hour  as  computed  by  the  AIM  model.  At  1000  local,  it  can  be  seen  that  there 
is  a  spike  in  the  maximum  value. 
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Table  4.  Surface  temperature  differences  (°C)  computed  using  different  solar  flux 
initialization  schemes  and  surface  conditions. 

Measured 

Surface  conditions 

Bare 

Snow 

Deciduous 

Coniferous 

Grass 

Maxdif 

1.0 

0.7 

0.6 

0.5 

0.4 

Mindif 

-0.2 

-0.2 

-0.2 

-0.1 

-0.1 

MODTRAN 

Avgdif 

0.4 

0.3 

0.2 

0.1 

0.2 

Maxdif 

0.0 

0.0 

0.0 

0.0 

0.0 

Mindif 

-1.8 

-1.3 

-1.1 

-0.7 

-0.7 

Shapiro 

Avgdif 

-0.8 

-0.5 

-0.5 

—0.3 

-0.3 

Maxdif 

0.4 

0.4 

0.2 

0.2 

0.2 

AIM 

Mindif 

-0.6 

-0.5 

-0.3 

—0.3 

-0.3 

Average 

Avgdif 

-0.1 

-0.1 

-0.1 

-0.1 

0.0 

Maxdif 

0.0 

0.0 

0.0 

0.0 

0.0 

AIM 

Mindif 

-4.0 

-2.9 

-2.2 

-1.3 

-1.5 

Maximum 

Avgdif 

-1.3 

—0.9 

-0.8 

-0.5 

-0.5 

Maxdif 

3.5 

2.7 

1.9 

1.2 

1.3 

AIM 

Mindif 

0.0 

0.0 

0.0 

0.0 

0.0 

Minimum 

Avgdif 

1.4 

1.1 

0.8 

0.5 

0.5 

Because  the  cloud  amount  at  this  time  is  0.8  low  cloud  and  0.2  middle  cloud, 
only  a  few  of  the  6400  pixels  have  high  flux  values.  Thus,  the  computed  average 
of  the  6400  values  tends  to  wash  out  the  peak  at  this  time.  The  measured  value  at 
1000  local  was  333  W/m2,  and  at  1009  it  was  only  96  W/m2.  Based  on  the  winds 
at  the  cloud  level,  in  9  minutes  the  clouds  could  have  advected  approximately  4 
km,  a  distance  significantly  larger  than  the  standard  gap  distance  between  clouds 
produced  in  the  AIM  domain.  This  is  a  limitation  of  nondistributed  flux  models 
or  point  measurements:  we  cannot  depict  the  spatial  and  temporal  variability  that 
can  occur  over  these  small  scales.  In  this  case,  it  just  happens  that  at  1000  local 
the  lack  of  cloud  cover  over  the  site  resulted  in  fairly  high  flux  values.  Another 
problem  that  comes  into  play  is,  how  do  models  stack  clouds  vertically?  Two 
approaches  are  random  overlap  and  maximum  overlap.  Stochastic  models  like 
CSSM  (the  model  used  in  this  work)  use  a  random  overlap  approach.  But  is 
nature  random?  If  there  are  dynamic  processes  involved,  it  is  likely  that  the 
clouds  will  be  correlated  from  level  to  level.  Based  on  the  AIM  maximum  and 
minimum  flux  initialization  schemes,  the  potential  range  of  temperatures  of  the 
AIM  domain  is  5°C  (Figure  12)  even  though  sky  conditions  are  cloudy.  Even  the 
range  of  temperatures  computed  using  MODTRAN,  average  AIM  values  over  the 
domain,  Shapiro’s  technique,  and  the  measured  is  on  the  order  of  several  degrees. 
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These  differences  reflect  the  difference  in  the  flux  values,  which  are  due  in  part 
to  how  clouds  are  handled  in  each  scheme.  During  the  period  from  1 000  to  1 500 
local,  all  schemes  predict  flux  values  that  are  greater  than  the  measured  values. 
Even  the  AIM  minimum  values  are  greater  than  the  measured  values.  This  is  not 
surprising,  considering  the  cloud  conditions.  The  cloud  types  were  mainly  cumu¬ 
lus  and  had  very  dark  bases  indicating  fairly  large  optical  depths.  As  pointed  out 
earlier,  information  on  the  vertical  extent  of  the  clouds  was  not  available.  Even  if 
this  information  were  available,  only  the  AIM  model  uses  cloud  thickness  infor¬ 
mation  in  the  calculation  of  the  solar  flux  values.  For  the  version  of  MODTRAN 
used  in  this  study  there  are  no  provisions  for  using  cloud  thickness  information. 
Since  the  coefficients  used  in  the  Shapiro  model  are  based  on  measurements, 
cloud  thickness  indirectly  plays  a  role  in  the  calculation  of  the  solar  flux  values. 


Yuma  1  -  26  Mar  93  (DOY  085) 
Total  Solar -Cloudy 


Figure  10.  Comparison  of  different  model  computed  flux  and  measured  flux 
for  cloudy  conditions  for  Yuma. 

The  impact  of  variations  in  the  cloud  conditions  is  reflected  in  the  summary 
statistics  presented  in  Table  5. 
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Figure  11.  AIM  maximum  fluxes  for  Yuma  cloudy  skies  day. 


Yuma  -  Bare  ground,  cloudy  skies  -  26  Mar  93  (DOY  085) 


Yuma  -  Cloudy  26  Mar  93  (DOY  085) 


Figure  12.  Surface  temperature  for  bare  ground  using  different  initialization 
schemes. 
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Table  5.  Summary  of  impact  of  different  solar  flux  initialization 
schemes  on  surface  temperature  (°C)  for  Yuma. 


Surface  conditions 

Measured 

Bare 

Grass 

Maxdif 

1.1 

0.5 

Mindif 

-2.0 

-0.6 

MODTRAN 

Avgdif 

-0.2 

-0.1 

Maxdif 

0.4 

0.4 

Mindif 

-5.0 

-1.6 

Shapiro 

Avgdif 

-15 

-0.4 

Maxdif 

3.9 

1.5 

AIM 

Mindif 

-3.2 

-1.7 

Average 

Avgdif 

-0.1 

0.0 

Maxdif 

0.2 

0.2 

AIM 

Mindif 

-6.9 

-2.4 

Maximum 

Avgdif 

-2.0 

-0.5 

Maxdif 

5.3 

1.9 

AIM 

Mindif 

-1.5 

-0.5 

Minimum 

Avgdif 

0.6 

0.2 

Grayling  II:  Cloudy  skies 

For  the  day  selected  as  representing  cloudy  conditions  during  the  Grayling  II 
experiment,  low  overcast  cloud  conditions  persisted  throughout  the  entire  day¬ 
light  period.  In  fact,  the  measured  direct  solar  flux  was  essentially  zero  for  all 
time  periods  except  around  local  noon.  Thus,  the  total  solar  loading  consisted 
only  of  the  diffuse  component.  MODTRAN  and  the  semi-empirical  model 
(Shapiro’s  model)  indicated  the  direct  component  was  zero  for  all  time  periods. 
A  plot  of  total  solar  flux  values  for  the  different  schemes  is  given  in  Figure  13. 

All  models  tend  to  overpredict  the  value  of  the  total  solar  flux  relative  to 
the  measured.  As  indicated,  the  total  solar  flux  is  essentially  the  diffuse  flux. 
Because  of  the  low  overcast  sky  conditions  we  have  no  information  on  cloud 
amounts  above  the  overcast.  If  upper  cloud  layers  (middle  and/or  high  clouds) 
were  present,  or  the  cloud  layer  was  optically  thicker  than  normal  for  this  type 
cloud,  the  optical  depth  would  increase,  and  the  corresponding  diffuse  solar 
component  may  decrease  if  the  total  optical  depth  is  greater  than  approximately 
3.  This  would  bring  the  model  values  more  in  line  with  the  measured  values. 

The  spike  in  the  measured  total  solar  flux  around  1300  local  is  due  to  a  direct 
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component  and  a  significant  increase  in  the  diffuse  component  (see  Figure  14). 
Recalling  Figure  4,  we  see  the  diffuse  component  can  increase  as  the  optical 
depth  increases  up  to  a  certain  point  and  then  decreases  as  the  optical  depth  con¬ 
tinues  to  increase.  In  addition,  the  diffuse  component  can  also  increase  due  to  an 
increase  in  the  scattering  from  the  sides  of  the  clouds.  Only  the  AIM  model  pre¬ 
dicted  a  direct  component  during  this  time  period.  A  histogram  of  the  6400  AIM 
values  and  50  measured  values  from  five  measurement  sites  for  a  period  of  5 
minutes  before  and  after  1300  local  are  presented  in  Figure  15.  The  measurement 
sites  are  within  a  circle  with  a  radius  on  the  order  of  1  km.  Even  over  this  small 
spatial  domain  for  this  relatively  short  period,  there  is  variability  in  the  flux 
values,  as  indicated  by  the  histogram  of  the  measured  values  in  Figure  15.  The 
impact  of  the  different  initialization  schemes  on  the  surface  temperature  for  bare 
ground  is  given  in  Figure  1 6.  As  anticipated,  the  greatest  range  of  temperature 
differences  occurs  around  solar  noon  when  SWOETHERM  is  initialized  using 
the  AIM  maximum  and  minimum  flux  values. 
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Grayling  II  -  Cloudy  skies  -  24  MAR  94  (DOY  083) 


Figure  13.  Total  solar  flux  values  for  overcast  low  cloud  conditions  at  Gray¬ 
ling  II  field  site. 

The  range  of  temperatures  for  the  other  surface  types  is  not  as  great  as  it  is 
for  bare  ground.  A  summary  of  the  results  for  the  different  initialization  schemes 
and  surface  material  types  is  given  in  Table  6. 
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Figure  14.  Comparison  of  direct  and  diffuse  measured  and  AIM-generated 
flux  values. 


Grayling  11-24  Mar  94  (DOY  083)  - 1 300  hr  -  Cloudy 
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Figure  15.  Histogram  of  total  solar  flux  for  6400  AIM  values  and  50  meas¬ 
ured  values. 
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Grayling  II  -  Bare  ground,  cloudy  skies  -  24  Mar  94  (DOY  083) 
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Figure  16.  Impact  of  different  flux  initialization  schemes  on  bare-ground 
surface  temperature  for  Grayling  II. 


Table  6.  Summary  of  surface  temperature  differences  (°C)  computed  using  different 
solar  flux  initialization  schemes  and  surface  conditions. 


Measured 


MODTRAN 


Shapiro 


AIM 

Average 


AIM 

Maximum 


AIM 

Minimum 
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In  almost  all  cases,  the  differences  between  the  temperature  computed  using 
the  measured  flux  values  and  either  MODTRAN,  Shapiro’s  technique,  or  the 
average  AIM  flux  values  are  negative,  implying  the  temperatures  computed  using 
the  model  flux  values  are  warmer  than  the  corresponding  temperatures  computed 
using  the  measured  flux  values.  In  general,  the  model  flux  values  are  greater  than 
the  measured  values.  This  highlights  the  importance  of  using  models  that  account 
for  the  differences  in  the  cloud  optical  depth  or  cloud  vertical  extent. 

Relying  on  surface  observations  to  characterize  the  cloud  conditions  has  a 
number  of  limitations.  First,  we  noted  cases  where  upper-level  cloud  layers  may 
have  been  present  but  were  obscured  from  the  view  of  the  observers.  Second, 
even  when  upper  levels  are  visible  through  gaps  in  the  lower  cloud  layers,  it  may 
not  be  possible  to  accurately  determine  the  fractional  cover  of  these  upper  layers. 
Finally,  cloud  optical  depth  or  cloud  vertical  extent  is  not  readily  available  from 
surface  observations.  Given  these  data  gaps,  other  methods  need  to  be  employed 
to  improve  the  accuracy  of  the  cloud  data.  Methods  might  include  improved,  or 
automated,  analysis  of  the  rawinsonde  data  to  determine  the  presence  of  cloud 
layers,  use  of  satellite  observations  to  determine  the  nature  of  upper  level  clouds, 
and/or  coupling  of  observer  information  with  outputs  from  numerical  weather 
predictions.  In  this  latter  case,  models  for  determining  positions  of  cloud  layers 
and  fractional  coverage  are  available,  such  as  the  Atmospheric  Sounding 
Program  (ASP)  (Passner  1998),  which  is  a  numerical  weather  post-processing 
algorithm  run  as  part  of  the  Army’s  Integrated  Meteorological  System  (IMETS). 
This  model  has  approximately  a  70%  skill  level  in  assessing  cloud  amount  and 
level  height  for  large  cloud  systems;  that  is,  clouds  not  involving  vertical  dynam¬ 
ics  like  cumulus  cloud  development.  Given  data  integration  of  this  sort,  it  would 
be  possible  to  gain  a  better  estimate  of  cloud  conditions.  Looking  at  the  problem 
from  another  perspective,  if  surface  solar  loading  information  were  available, 
model  outputs  could  be  used  to  infer  cloud  conditions  above  an  obscuring  lower 
cloud  level. 

Partly  cloudy  sky  flux  analysis 

As  indicated  in  the  introduction,  distributed  models  potentially  require 
distributed  solar  flux  values,  especially  for  partly  cloudy  conditions.  In  partly 
cloudy  conditions,  low  clouds  cast  hard  shadows  and  it  is  possible  to  transition 
from  full  solar  loading  conditions  to  diffuse  solar  loading  over  small  spatial 
extents  and  short  periods  of  time.  The  solar  loading  at  any  given  point  is  modu¬ 
lated  by  the  advection  of  clouds  and,  to  some  extent,  the  evolution  of  the  cloud 
geometry.  In  dry  environments  such  as  Yuma,  there  can  be  dramatic  changes 
in  surface  temperature  over  short  periods  of  time  or  over  small  spatial  extents. 
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Figure  1 7  clearly  illustrates  the  impact  of  temporal  variations  in  the  measured 
total  solar  flux  on  the  measured  surface  temperature. 

Yuma  Site  D  -27  March 


Surface  Temperature 


Figure  17.  Solar  flux  and  surface  temperature  for  site  D,  27  March,  Yuma. 

In  a  12-min  period  between  1337  (13.6)  and  1349  (13.8)  the  surface  tempera¬ 
ture  dropped  3.8°C  and  the  solar  flux  changed  from  approximately  900  W/m2  to 
approximately  240  W/m2.  These  fluctuations  in  the  total  solar  flux  are  associated 
with  partly  cloudy  sky  conditions.  In  dry,  desert-like  climates,  the  surface  tem¬ 
perature  responds  rapidly  to  changes  in  solar  loading.  Figure  1 8  depicts  the  spa¬ 
tial  and  temporal  variability  in  the  surface  temperature  and  the  total  solar  flux  for 
sites  C  and  D.  Site  C  was  located  in  an  area  characterized  by  a  hard-packed,  dark- 
colored  surface.  Site  D  was  located  in  the  same  area,  but  it  was  approximately 
1 00  m  south  of  site  D.  The  ratio  of  the  surface  temperature  at  site  C  to  that  at  site 
D  is  less  than  1 ,  for  the  3-hr  period  depicted  in  Figure  1 8.  It  is  believed  this  is 
due  to  differences  between  the  two  sites  (slope,  vegetation  cover,  soil  type,  etc.). 
The  ratio  of  the  total  solar  flux  for  the  two  sites  is  also  depicted  in  this  figure.  At 
times  this  ratio  has  values  as  high  as  2,  indicating  the  total  solar  flux  at  one  site  is 
twice  that  at  the  other  site,  even  though  the  two  sites  are  only  about  100  m  apart. 
There  are  also  periods  when  the  ratio  is  1 .  An  inspection  of  the  data  reveals  that 
most  of  these  periods  correspond  to  clear  conditions.  For  the  most  part,  the  sur¬ 
face  temperature  ratio  is  in  phase  with  the  solar  flux  ratio  (there  is  a  small  time 
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lag).  Two  exceptions  are  around  1330  (13.5)  and  after  1530  (15.5).  During  the 
first  period  the  solar  fluxes  at  both  sites  are  on  the  order  of  250  W/m2,  indicating 
there  is  little  or  no  direct  component  of  solar  flux.  During  the  second  period, 
conditions  are  clear  and  the  total  solar  flux  consists  of  both  the  direct  and  diffuse 
component. 


Figure  18.  Surface  temperature  and  solar  flux  ratio  for  sites  C  and  D. 

While  it  cannot  be  verified,  it  is  believed  that  the  variations  in  the  surface 
temperature  ratio  during  these  periods  reflect,  in  part,  variations  in  the  other 
components  of  the  energy  balance. 

Grayling  I:  Partly  cloudy  skies 

During  the  daylight  hours,  the  cloud  cover  varied  from  nine-tenths  coverage 
during  the  early  morning  to  one-tenth  coverage  in  the  late  afternoon. 

The  distribution  of  the  total  solar  flux  for  both  the  measured  values  from  site 
El  and  E2  and  the  AIM-calculated  values  is  presented  in  Figure  19. 

The  measured  values  are  the  1-min  values  for  the  daylight  period  for  loca¬ 
tions  El  and  E2.  The  AIM  values  are  the  values  calculated  for  an  8-  by  8-km  area 
based  on  the  cloud  amounts  reported  on  the  hour.  Thus,  the  AIM  values  are  a 
snapshot  over  an  8-  by  8-km  aerial  extent  on  the  hour,  while  the  measured  values 
are  the  minute-by-minute  values  for  a  single  location.  AIM  has  a  higher  percent¬ 
age  of  the  flux  values  in  the  lower  bin  ranges.  Even  the  measured  values  from 
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two  sites,  separated  by  approximately  100  m,  show  significant  differences  in  the 
percent  occurrence  of  flux  values.  On  minute-by-minute  bases,  the  shadowing  at 
the  two  sites  varies  significantly. 
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Figure  19.  Distribution  of  total  solar  flux  for  sites  El  and  E2  for  Grayling  I 
and  AIM-calculated  values. 

From  Table  7  it  is  evident  that  the  hourly  reported  cloud  cover  changed 
during  the  day.  In  general,  there  was  a  clearing  trend.  Potentially,  the  AIM 
distribution  of  flux  values  may  more  closely  resemble  the  measured  values  if  the 
model  was  run  more  frequently.  Unfortunately,  it  is  not  possible  to  do  this  since 
cloud  observations  were  not  taken  more  frequently  than  once  per  hour.  These 
flux  values  along  with  the  flux  values  computed  from  Shapiro’s  technique  and 
MODTRAN  have  been  used  to  ascertain  the  response  of  the  surface  temperature. 
Figure  20  presents  the  results  for  bare  ground,  and  Table  8  summarizes  the 
results  for  all  surfaces.  As  anticipated,  the  surface  temperatures  calculated  using 
the  AIM  maximum  and  minimum  values  of  solar  flux  for  each  hour  bracket  the 
other  values.  As  pointed  out  earlier,  we  would  not  anticipate  experiencing  this 
range  of  temperatures  because  a  single  spatial  location  would  not  experience 
either  the  minimum  or  the  maximum  solar  value  for  the  entire  day.  What  it  does 
indicate  is  that  we  can  experience  a  range  of  temperatures  at  a  single  location 
due  to  variations  in  solar  loading.  In  fact,  the  range  of  temperature  differences 
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exceeds  the  stated  accuracy  of  the  thermal  model.  In  modeling  synthetic  scenes, 
it  is  important  either  to  model  accurately  the  solar  flux  values  or  to  state  the 
results  in  terms  of  a  potential  temperature  range  rather  than  a  single  value. 


Table  7.  Cloud  amounts  used  in  solar  flux  models 
for  partly  cloudy  conditions  for  Grayling  1. 

Cloud  amount  (tenths) 

Time 

Low 

Middle 

High 

6 

0.0 

0.0 

0.0 

7 

0.0 

0.6 

0.0 

8 

0.7 

0.1 

0.0 

9 

0.9 

0.0 

0.0 

10 

0.9 

0.0 

0.0 

11 

0.8 

0.0 

0.0 

■a 

0.6 

0.0 

0.0 

■9 

0.6 

0.0 

0.0 

14 

0.3 

0.0 

0.0 

15 

0.3 

0.0 

0.0 

16 

0.1 

0.0 

17 

0.1 

1 

0.1 

18 

0.1 

MiSEKl 

0.1 

19 

0.1 

0.0 

B : : 

20 

0.1 

0.0 

1 

Yuma:  Partly  cloudy  skies 

The  variations  in  the  cloud  cover  are  given  in  Table  9.  Maximum  cloud  cover 
occurs  around  solar  noon. 

In  Figure  2 1 ,  the  distribution  of  the  AIM  total  solar  flux  at  1 300  local  is 
depicted.  As  indicated  in  Table  9,  the  total  cloud  amount  at  1 300  local  is  70%. 
Flux  values  greater  than  1000  W/m2  cover  only  7%  of  the  8-  by  8-km  area  at 
1300  local,  but  51%  of  the  area  has  flux  values  between  500  and  1000  W/m2, 
with  the  majority  of  these  values  in  the  500  to  600  W/m2  range.  Flux  values  in 
the  400  to  500  W/m2  range  cover  another  22%  of  the  area.  In  this  dry  climate,  the 
surface  temperature  will  respond  rapidly  to  the  wide  range  of  flux  values.  As  in¬ 
dicated  earlier,  the  measured  surface  temperature  changed  approximately  4.0°C 
in  just  12  minutes. 
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GRAYLING  I  -  Bare  ground,  partly  cloudy  skies  -  24  Oct  92  (DOY  298) 
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Figure  20.  Surface  temperature  for  bare  ground  resulting  from  different 
solar  flux  initialization  schemes. 


Table  8.  Summary  of  surface  temperature  differences  (°C)  computed  using  different 
solar  flux  initialization  schemes  and  surface  conditions. 

Measured 

Surface  conditions 

Bare 

Snow 

Deciduous 

Coniferous 

Grass 

Maxdif 

MSM 

1.9 

1.3 

Kwm 

Mindif 

0.0 

0.0 

MODTRAN 

Avgdif 

ISM 

1.0 

0.9 

0.6 

0.5 

1.4 

mum 

1.0 

0.7 

0.5 

-0.2 

“0.1 

“0.1 

Shapiro 

:-:Wm 

0.2 

0.4 

0.3 

0.1 

4.6 

3.2 

2.5 

msm 

Mm 

AIM 

0.0 

0.0 

0.0 

■ 

Average 

2.3 

1.8 

1.0 

0.6 

0.7 

■ 

1.2 

0.9 

HBl 

AIM 

1 

-0.9 

“0.4 

Maximum 

0.0 

0.0 

0.0 

0.0 

7.0 

5.3 

3.7 

2.1 

WEm 

AIM 

0.0 

0.0 

0.0 

0.0 

wBM 

Minimum 

3.3 

2.7 

1.5 

0.8 

i.i 
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Table  9.  Cloud  amounts  used  in  solar  flux  models. 

Cloud  amount  (tenths) 

Time 

Low 

Middle 

High 

7 

0.2 

0.0 

0.0 

8 

0.5 

0.0 

0.0 

9 

0.0 

0.0 

0.0 

10 

0.0 

0.0 

0.0 

11 

0.1 

0.0 

0.0 

12 

0.1 

0.0 

0.0 

13 

0.7 

0.0 

0.0 

14 

0.6 

0.0 

0.0 

15 

0.7 

0.0 

0.0 

16 

0.4 

0.0 

0.2 

17 

0.3 

0.0 

0.3 

18 

0.1 

0.1 

0.1 
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Figure  21.  Distribution  of  total  solar  flux  at  1300  local  as  calculated  by  AIM 
model. 
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Figure  22.  Mean  AIM-calculated  total  solar  flux  ±  one  standard  deviation. 

The  standard  deviation  associated  with  the  temporal  variation  of  the  meas¬ 
ured  flux  and  the  spatial  variation  of  the  AIM  flux  values  is  given  in  Figure  22. 
For  the  AIM  total  solar  flux,  this  standard  deviation  is  based  on  the  6400  values 
in  the  8-  by  8-km  area,  and  for  the  measured  value  it  is  the  deviation  associated 
with  the  60  previous  1-min  values.  The  standard  deviation  at  0900  and  1000  local 
is  0  for  the  AIM  values  since  there  is  no  cloud  cover  at  these  times.  The  largest 
standard  deviation  for  both  the  measured  data  and  the  AIM  calculated  values 
occurs  at  1400  local.  The  standard  deviation  associated  with  the  temporal  vari¬ 
ation  (measured  data)  is  approximately  313  W/m2;  the  standard  deviation  asso¬ 
ciated  with  the  spatial  variability  (AIM  calculated  values)  is  approximately  263 
W/m2.  If  we  assume  no  evolution  in  the  cloud  geometry,  then  there  should  be  a 
relationship  between  the  spatial  and  temporal  variation  of  the  solar  flux.  At  this 
point  in  time  we  do  not  have  sufficient  information  to  establish  this  relationship, 
and  the  information  in  Figure  22  is  presented  to  emphasize  the  spatial  and  tem¬ 
poral  variability  of  the  solar  flux  for  partly  cloudy  conditions  and  the  resulting 
impact  on  the  surface  temperature.  The  spatial  and  temporal  variations  in  the 
surface  temperature  represent  clutter  to  an  infrared  sensor  system  and  may  reduce 
the  performance  of  that  system.  The  impact  of  the  different  initialization  schemes 
on  a  vegetated  surface  is  depicted  in  Figure  23.  The  resulting  surface  temperature 
from  all  schemes  tends  to  follow  the  surface  temperature  as  predicted  by  using 
the  measured  fluxes.  In  addition,  the  range  of  vegetation  temperature  values  as 
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predicted  from  using  the  maximum  and  minimum  is  not  as  great  as  for  bare 
ground  (Table  10).  The  average  difference  between  the  surface  temperature 
computed  using  the  measured  flux  and  the  surface  temperature  using  the  AIM 
maximum  flux  is  -3.1°C,  indicating  AIM  maximum  fluxes  result  in  warmer 
surface  temperatures.  The  average  difference  using  the  AIM  minimum  flux  is 
4.6°C.  Thus,  the  range  of  surface  temperatures  based  on  the  average  of  the 
temperatures  computed  using  the  maximum  and  minimum  AIM  fluxes  is  7.7°C. 


Yuma  MVEG  -  Partly  cloudy  skies  - 14  Apr  93  (DOY 104) 


Figure  23.  Surface  temperature  for  medium  vegetation  (grass)  for  different 
solar  flux  initialization  schemes. 

During  the  1300  to  1400  local  time  period,  the  measured  surface  temperature 
difference  between  sites  C  and  D  was  as  high  as  2.0°C.  These  sites  are  separated 
by  less  than  100  m.  During  the  night,  the  surface  temperature  difference  is  on  the 
order  of  0.2  to  0.3°C.  The  differences  during  the  1300  to  1400  local  time  period 
are  believed  to  be  due  to  differential  solar  loading  resulting  from  the  partly 
cloudy  conditions. 

Grayling  II:  Partly  cloudy  skies 

Before  noon  local  time,  sky  conditions  ranged  from  partly  cloudy  to  over¬ 
cast.  During  the  afternoon,  skies  cleared.  The  changes  in  the  cloud  conditions 
and  their  impact  on  the  total  solar  flux  and  the  surface  temperature  are  given  in 
Table  1 1  and  Figure  24,  respectively. 
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Table  10.  Summary  of  surface  temperature  differences  (°C)  computed 
using  different  solar  flux  initialization  schemes  and  surface  conditions. 

Surface  conditions 

Measured 

Bare 

Grass 

5.8 

2.0 

■ 

-1.4 

-0.6 

MODTRAN 

| 

0.9 

0.3 

■ 

2.0 

0.7 

1 

-5.2 

-1.9 

Shapiro 

| 

-0.3 

-0.1 

■ 

6.5 

2.2 

AIM 

1 

-0.5 

-0.4 

Average 

1.8 

0.6 

0.4 

0.2 

AIM 

1 

-10.1 

-3.5 

Maximum 

-1.3 

-0.4 

14.7 

4.6 

AIM 

Mindif 

0.0 

0.0 

Minimum 

Avgdif 

4.6 

1.4 

The  total  solar  flux  during  the  morning  varies  by  a  factor  of  three  over 
relatively  short  time  periods.  Not  only  does  the  surface  temperature  follow  the 
general  trend  of  the  solar  flux  but  it  also  reflects  some  of  the  high-frequency 
changes  in  the  solar  flux.  Changes  in  surface  temperature  occur  on  the  order  of 
several  degrees.  At  1000  local,  the  sky  conditions  were  basically  overcast  with 
low  cloud.  Under  these  conditions,  AIM  calculated  total  solar  flux  values  in  a 
narrow  distribution  with  values  that  are  typical  of  diffuse  solar  radiation  (Figure 
25).  On  the  other  hand,  the  measured  values  are  considerably  higher  and  spread 
over  a  greater  range.  The  AIM  fluxes  are  spatially  distributed  fluxes  based  on  the 
cloud  cover  at  a  single  time,  and  the  measured  fluxes  are  temporally  distributed 
over  the  previous  60  minutes  (60  observations).  In  this  case,  it  appears  the  spatial 
distribution  cannot  be  equated  to  the  temporal  distribution. 

The  AIM-calculated  total  solar  flux  in  Figure  26  is  the  average  for  the  6400 
values  for  the  8-  by  8-km  area  based  on  the  cloud  amount  on  the  hour.  In  general, 
the  AIM  total  solar  flux  is  less  during  cloud  periods,  but  rapidly  converges  to  the 
other  values  when  the  skies  clear  after  1500  local.  The  spatial  distribution  of  the 
clouds  is  based  on  CSSM.  CSSM  does  have  a  capability  to  produce  the  temporal 
variation  of  clouds  for  short  periods  of  time  (on  the  order  of  1 5  min)  but  running 
AIM  for  minute-by-minute  cloud  distributions  would  be  computation-prohibitive 


56 


ERDC/CRREL  TR-03-13 


at  this  time.  Accurate  high-resolution  modeling  requires  both  the  spatial  and  tem¬ 
poral  distribution  of  clouds  and  the  associated  solar  fluxes.  The  calculated  sur¬ 
face  temperature  based  on  the  different  initialization  schemes  is  given  in  Figure 
27.  In  the  afternoon  when  skies  clear,  the  surface  temperatures  computed  using 
the  different  initialization  schemes  are  almost  identical.  Earlier  in  the  day  when 
it  is  cloudy  this  is  not  the  case.  In  fact,  the  temperature  difference  is  as  large  as 
10°C  for  the  initialization  schemes  based  on  the  maximum  and  minimum  AIM- 
calculated  flux.  In  general,  the  values  given  in  Table  12  for  the  different  initiali¬ 
zation  schemes  and  surface  conditions  are  greater  than  the  corresponding  values 
for  overcast  and  clear-sky  conditions.  Partly  cloudy  skies  will  result  in  large 
spatial  and  temporal  variations  in  solar  loading. 


Table  11.  Cloud  amounts  used  in  solar  flux  models. 

Cloud  amount  (tenths) 

Time 

Low 

Middle 

High 

6 

0.1 

0.2 

■ai 

7 

0.8 

0.0 

8 

1.0 

0.0 

lEl 

9 

1.0 

0.0 

0.0 

10 

1.0 

0.0 

0.0 

11 

0.9 

0.0 

0.1 

mm 

0.6 

0.0 

mm 

0.4 

0.0 

PI 

14 

0.2 

0.0 

15 

0.1 

0.0 

0.1 

16 

mmrma 

0.0 

17 

0.0 

W3X 

18 

Ea 

0.0 

0.0 

19 

0.0 

0.0 

0.0 

20 

0.0 

0.0 

0.0 
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Figure  24.  Total  solar  flux  and  surface  temperature  measurements  at  site 
E3  for  partly  cloudy  conditions. 
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Figure  25.  Distribution  of  measured  total  solar  flux  for  period  of  60  minutes 
for  Grayling  II  and  AIM-calculated  values  for  8-  by  8-km  area  at  1000  local. 
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Grayling  II  -  10  Apr  94  (DOY  100) 


Figure  26.  Total  solar  flux  measured  and  calculated  for  partly  cloudy  condi¬ 
tions  for  Grayling  II. 


Grayling  II  -  Bare  ground,  partly  cloudy  skies  - 10  Apr 94  (DOY  100) 


Figure  27.  Calculated  surface  temperature  based  on  different  initialization 
schemes. 
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5  DISCUSSION  AND  CONCLUSION 

Distributed  energy  budget  models  require  distributed  initialization  condi¬ 
tions.  The  optimal  response  of  snowmelt  using  the  one-dimensional  SNTHERM 
(the  basic  model  used  in  SWOETHERM)  model  in  a  distributed  mode  for  mete¬ 
orological  cell  sizes  ranging  from  9  km  to  72  km  for  the  Sava  River  Basin  was 
obtained  using  a  cell  size  of  9  km.  IMETS  is  being  designed  to  provide  meteoro¬ 
logical  information  over  spatial  extents  as  large  as  500-  by  500-km  with  a  cell 
resolution  of  approximately  10  km  (higher  resolutions  are  being  planned  for  the 
future).  The  response  of  the  surface  temperature  for  cloudy  and  partly  cloudy 
conditions  can  be  on  the  order  of  several  degrees  for  temporal  and  spatial  varia¬ 
tions  in  the  solar  flux  on  the  order  of  tens  of  minutes  and  several  hundred  meters. 
Of  the  solar  flux  models  investigated,  only  AIM  provides  the  spatial  distribution 
of  the  solar  flux.  This  model  also  has  the  potential  to  provide  the  temporal  varia¬ 
tions  of  the  solar  flux  if  it  is  coupled  with  a  cloud  model  that  generates  both  the 
spatial  and  temporal  variation  in  cloud  cover.  CSSM  has  this  capability.  Run 
times  of  models  like  AIM  that  can  provide  spatial,  temporal,  and  spectral  infor¬ 
mation  are  an  issue.  We  have  investigated  a  number  of  techniques  to  decrease  the 
run  time  of  AIM,  and  additional  acceleration  techniques  will  be  explored  in  the 
future.  For  clear-sky  conditions  and  small  spatial  extents,  semi-empirical  models 
like  Shapiro’s  model  provide  sufficiently  accurate  flux  information.  The  issue 
becomes  more  complicated  when  we  consider  overcast  or  partly  cloudy  skies. 
Even  for  overcast  conditions,  we  observed  both  spatial  and  temporal  variations 
in  the  measured  solar  fluxes  that  resulted  in  surface  temperature  variations  of 
several  degrees.  These  variations  in  the  solar  flux  have  long  been  attributed  to 
variations  in  the  cloud  optical  properties  (Welch  and  Wielicki  1984,  Li  et  al. 
1994).  We  have  shown  that  as  the  aerosol  optical  depth  increases,  the  direct 
component  decreases  in  a  nonlinear  fashion  and  the  diffuse  component  first 
increases  and  then  decreases.  We  would  expect  a  similar  behavior  for  variations 
in  cloud  optical  depths.  This  implies  that  in  addition  to  knowing  the  horizontal 
distribution  of  clouds  it  is  also  necessary  to  know  typical  cloud  optical  depths. 

As  indicated,  the  aerosol  optical  depth  can  impact  the  surface  solar  flux. 
Frequently  the  visibility  is  used  as  a  surrogate  for  the  aerosol  optical  depth.  What 
is  really  needed  is  the  variation  of  the  optical  depth  with  height.  For  clear  and 
partly  cloudy  conditions,  the  surface  predicted  temperature  based  on  the  solar 
flux  predicted  by  Shapiro’s  model  follows  closely  the  values  predicted  using  the 
measured  fluxes.  The  predictions  are  not  as  accurate  for  cloudy  conditions,  and 
thick  haze  conditions  were  not  studied.  While  Shapiro’s  technique  does  allow 
for  different  cloud  types,  it  does  not  allow  for  variations  in  the  cloud  thickness 
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(variations  in  the  cloud  optical  depth).  In  addition,  the  predictions  based  on 
Shapiro’s  model  are  for  a  single  location  at  a  single  time.  This  model  does  not 
provide  for  spatial  variations  in  the  solar  flux. 

Another  potential  issue  with  respect  to  Shapiro’s  model  is  the  relationship 
between  observer-collected  cloud  cover  information  and  the  model-computed 
surface  solar  flux.  Due  to  automated  analysis  techniques  being  developed  for 
satellite  or  lidar  retrieval  of  cloud  data,  it  may  be  that  the  availability  of  human- 
observed  cloud  condition  information  will  be  at  a  premium  in  the  future.  In  this 
case,  questions  need  to  be  raised  as  to  the  relationship  between  Shapiro’s  coef¬ 
ficients  and  the  actual  cloud  state.  Shapiro’s  model  appears  to  have  a  good  skill 
level  in  inferring  solar  flux  given  observed  cloud  conditions.  Hence,  whatever 
correlations  are  present  between  observed  conditions  and  observed  fluxes  are 
accurately  being  preserved  through  his  model.  Thus,  the  effects  of  clouds  that 
are  unobservable  are  being  convolved  into  the  model  responses,  and  correlations 
between  these  unseen  upper  cloud  layers  are  included  in  the  model  coefficients 
applied  strictly  to  the  lower  layer  effects.  These  effects  might  explain  the  model’s 
(Shapiro’s)  superior  performance  in  handling  some  of  the  cloudy  cases.  How¬ 
ever,  if  the  method  of  measuring  cloud  conditions  changes,  the  Shapiro  model 
would  need  to  be  altered,  since  the  same  correlations  might  not  exist  when  differ¬ 
ent  data  acquisition  methods  are  used.  One  area  of  active  exploration  might  be  to 
inquire,  using  a  3-D  model,  what  cloud  optical  conditions  would  be  required  to 
produce  specific  outputs  from  the  Shapiro  model.  These  considerations  could 
lead  to  extensions  of  the  Shapiro  model  under  various  cloud  conditions. 

Spatial  and  temporal  variations  in  the  state  of  the  ground  can  negatively 
impact  weapon  system  performance.  Sufficient  spatial  variation  in  the  surface 
temperature  can  confound  infrared  target  detection.  And,  as  noted  in  the  data 
analysis,  often  measurements  taken  less  than  a  kilometer  apart  showed  significant 
variations  in  the  sensed  incident  fluxes.  Variations  in  snowmelt  conditions  due 
to  variations  in  solar  loading  can  negatively  impact  active  radar  systems  as  well. 
Techniques  are  being  developed  to  segment  terrain  in  order  to  run  one-dimen¬ 
sional  energy  budget  models  over  aerial  extents  for  tens  of  kilometers  to  portray 
the  battlespace  state  of  the  ground  and  the  impact  of  the  environment  on  military 
systems.  Accurate  depiction  of  the  state  of  the  ground  for  synthetic  scene  simu¬ 
lations  depends,  in  part,  on  the  terrain  segmentation  techniques  employed  and  the 
accuracy  of  the  distributed  meteorological  information  used  as  boundary  condi¬ 
tions  for  the  energy  budget  models  that  calculate  the  state  of  the  ground. 
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